Line data Source code
1 : /*
2 : * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__FUNCTION_SPACES__INTEGRATE__
34 : #define __H__UG__LIB_DISC__FUNCTION_SPACES__INTEGRATE__
35 :
36 : #include <cmath>
37 :
38 : #include <boost/function.hpp>
39 :
40 : #include "common/common.h"
41 :
42 : #include "lib_grid/tools/subset_group.h"
43 :
44 : #include "lib_disc/common/function_group.h"
45 : #include "lib_disc/common/groups_util.h"
46 : #include "lib_disc/domain_util.h" // for CollectCornerCoordinates
47 : #include "lib_disc/quadrature/quadrature_provider.h"
48 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
49 : #include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
50 : #include "lib_disc/spatial_disc/user_data/user_data.h"
51 : #include "lib_disc/spatial_disc/user_data/const_user_data.h"
52 : #include "lib_disc/reference_element/reference_mapping_provider.h"
53 :
54 : #ifdef UG_FOR_LUA
55 : #include "bindings/lua/lua_user_data.h"
56 : #endif
57 :
58 : namespace ug{
59 :
60 :
61 : ////////////////////////////////////////////////////////////////////////////////
62 : // Object encapsulating (scalar) GridFunction related data
63 : ////////////////////////////////////////////////////////////////////////////////
64 : template <typename TGridFunction>
65 : class ScalarGridFunctionData
66 : {
67 : public:
68 : typedef typename TGridFunction::domain_type domain_type;
69 :
70 0 : ScalarGridFunctionData(TGridFunction& gridFct, size_t cmp)
71 0 : : m_gridFct(gridFct), m_fct(cmp),
72 0 : m_id(m_gridFct.local_finite_element_id(m_fct))
73 : {}
74 :
75 : TGridFunction& grid_function()
76 0 : {return m_gridFct;}
77 :
78 : const TGridFunction& grid_function() const
79 : {return m_gridFct;}
80 :
81 : size_t fct()
82 0 : {return m_fct;}
83 :
84 : const LFEID &id() const
85 0 : {return m_id;}
86 :
87 : //! returns true, iff scalar function is defined in subset si
88 : bool is_def_in_subset(int si) const
89 0 : { return m_gridFct.is_def_in_subset(m_fct, si); }
90 :
91 : /// returns domain (forward)
92 0 : SmartPtr<domain_type> domain() {return m_gridFct.domain();}
93 :
94 : /// returns const domain (forward)
95 : ConstSmartPtr<domain_type> domain() const {return m_gridFct.domain();}
96 :
97 : template <typename TElem>
98 : size_t dof_indices(TElem* elem, std::vector<DoFIndex>& ind, bool bHang = false, bool bClear = true) const
99 0 : {return m_gridFct.dof_indices(elem, m_fct, ind, bHang, bClear);}
100 :
101 : private:
102 : /// grid function
103 : TGridFunction& m_gridFct;
104 :
105 : /// component of function
106 : size_t m_fct;
107 :
108 : /// local finite element id
109 : LFEID m_id;
110 : };
111 :
112 : ////////////////////////////////////////////////////////////////////////////////
113 : // Integrand Interface
114 : ////////////////////////////////////////////////////////////////////////////////
115 :
116 :
117 : /// Abstract integrand interface
118 : /*! An integrand is a short-living (temporary) object that is generated/used in various integration functions.*/
119 : template <typename TData, int TWorldDim>
120 : class IIntegrand
121 : {
122 : public:
123 : /// world dimension
124 : static const int worldDim = TWorldDim;
125 :
126 : /// data type
127 : typedef TData data_type;
128 :
129 : /// returns the values of the integrand for a bunch of ips
130 : /**
131 : *
132 : * @param vValue[out] the value of the integrand at the ips
133 : * @param vGlobIP[in] global integration point positions
134 : * @param pElem[in] the element to integrate
135 : * @param vCornerCoords[in] corner coordinates of the element
136 : * @param vLocIP[in] local integration point positions
137 : * @param vJT[in] jacobian transposed at integration point
138 : * @param numIP[in] number of integration points
139 : */
140 : /// \{
141 : virtual void values(TData vValue[],
142 : const MathVector<worldDim> vGlobIP[],
143 : GridObject* pElem,
144 : const MathVector<worldDim> vCornerCoords[],
145 : const MathVector<1> vLocIP[],
146 : const MathMatrix<1, worldDim> vJT[],
147 : const size_t numIP) = 0;
148 : virtual void values(TData vValue[],
149 : const MathVector<worldDim> vGlobIP[],
150 : GridObject* pElem,
151 : const MathVector<worldDim> vCornerCoords[],
152 : const MathVector<2> vLocIP[],
153 : const MathMatrix<2, worldDim> vJT[],
154 : const size_t numIP) = 0;
155 : virtual void values(TData vValue[],
156 : const MathVector<worldDim> vGlobIP[],
157 : GridObject* pElem,
158 : const MathVector<worldDim> vCornerCoords[],
159 : const MathVector<3> vLocIP[],
160 : const MathMatrix<3, worldDim> vJT[],
161 : const size_t numIP) = 0;
162 : /// \}
163 :
164 : virtual ~IIntegrand() {}
165 :
166 :
167 : /// sets the subset
168 0 : virtual void set_subset(int si) {m_si = si;}
169 :
170 : /// returns the subset
171 0 : inline int subset() const {return m_si;}
172 :
173 : protected:
174 : /// subset
175 : int m_si;
176 : };
177 :
178 : /// Abstract integrand interface (using CRTP)
179 : template <typename TData, int TWorldDim, typename TImpl>
180 0 : class StdIntegrand : public IIntegrand<TData, TWorldDim>
181 : {
182 : public:
183 : /// world dimension
184 : static const int worldDim = TWorldDim;
185 :
186 : /// data type
187 : typedef TData data_type;
188 :
189 : /// \copydoc IIntegrand::values
190 : /// \{
191 0 : virtual void values(TData vValue[],
192 : const MathVector<worldDim> vGlobIP[],
193 : GridObject* pElem,
194 : const MathVector<worldDim> vCornerCoords[],
195 : const MathVector<1> vLocIP[],
196 : const MathMatrix<1, worldDim> vJT[],
197 : const size_t numIP)
198 : {
199 0 : getImpl().template evaluate<1>(vValue,vGlobIP,pElem,vCornerCoords,vLocIP,vJT,numIP);
200 0 : }
201 0 : virtual void values(TData vValue[],
202 : const MathVector<worldDim> vGlobIP[],
203 : GridObject* pElem,
204 : const MathVector<worldDim> vCornerCoords[],
205 : const MathVector<2> vLocIP[],
206 : const MathMatrix<2, worldDim> vJT[],
207 : const size_t numIP)
208 : {
209 0 : getImpl().template evaluate<2>(vValue,vGlobIP,pElem,vCornerCoords,vLocIP,vJT,numIP);
210 0 : }
211 0 : virtual void values(TData vValue[],
212 : const MathVector<worldDim> vGlobIP[],
213 : GridObject* pElem,
214 : const MathVector<worldDim> vCornerCoords[],
215 : const MathVector<3> vLocIP[],
216 : const MathMatrix<3, worldDim> vJT[],
217 : const size_t numIP)
218 : {
219 0 : getImpl().template evaluate<3>(vValue,vGlobIP,pElem,vCornerCoords,vLocIP,vJT,numIP);
220 0 : }
221 : /// \}
222 :
223 : protected:
224 : /// access to implementation
225 : TImpl& getImpl() {return static_cast<TImpl&>(*this);}
226 :
227 : /// const access to implementation
228 : const TImpl& getImpl() const {return static_cast<const TImpl&>(*this);}
229 : };
230 :
231 : ////////////////////////////////////////////////////////////////////////////////
232 : ////////////////////////////////////////////////////////////////////////////////
233 : // Generic Volume Integration Routine
234 : ////////////////////////////////////////////////////////////////////////////////
235 : ////////////////////////////////////////////////////////////////////////////////
236 :
237 : /// integrates on the whole domain
238 : /**
239 : * This function integrates an arbitrary integrand over the whole domain.
240 : * Note:
241 : * - only grid elements of the same dimension as the world dimension of the
242 : * domain are integrated. Thus, no manifolds.
243 : * - The implementation is using virtual functions. Thus, there is a small
244 : * performance drawback compared to hard coding everything, but we gain
245 : * flexibility. In addition all virtual calls compute for the whole set of
246 : * integration points to avoid many virtual calls, i.e. only one virtual
247 : * call for all integration points is needed.
248 : *
249 : * \param[in] iterBegin iterator to first geometric object to integrate
250 : * \param[in] iterBegin iterator to last geometric object to integrate
251 : * \param[in] integrand Integrand
252 : * \param[in] quadOrder order of quadrature rule
253 : * \param[in] quadType
254 : * \param[in] paaElemContribs (optional). If != NULL, the method will store
255 : * the contribution of each element in the
256 : * associated attachment entry.
257 : * \returns value of the integral
258 : */
259 : template <int WorldDim, int dim, typename TConstIterator>
260 0 : number Integrate(TConstIterator iterBegin,
261 : TConstIterator iterEnd,
262 : typename domain_traits<WorldDim>::position_accessor_type& aaPos,
263 : IIntegrand<number, WorldDim>& integrand,
264 : int quadOrder, std::string quadType,
265 : Grid::AttachmentAccessor<
266 : typename domain_traits<dim>::grid_base_object, ANumber>
267 : *paaElemContribs = NULL
268 : )
269 : {
270 : PROFILE_FUNC();
271 :
272 : // reset the result
273 : number integral = 0.0;
274 :
275 : // note: this iterator is for the base elements, e.g. Face and not
276 : // for the special type, e.g. Triangle, Quadrilateral
277 : TConstIterator iter = iterBegin;
278 :
279 : // this is the base element type (e.g. Face). This is the type when the
280 : // iterators above are dereferenciated.
281 : typedef typename domain_traits<dim>::grid_base_object grid_base_object;
282 :
283 : // get quad type
284 0 : if(quadType.empty()) quadType = "best";
285 0 : QuadType type = GetQuadratureType(quadType);
286 :
287 : // accessing without dereferencing a pointer first is simpler...
288 : Grid::AttachmentAccessor<grid_base_object, ANumber> aaElemContribs;
289 0 : if(paaElemContribs)
290 0 : aaElemContribs = *paaElemContribs;
291 :
292 : // We'll reuse containers to avoid reallocations
293 : std::vector<MathVector<WorldDim> > vCorner;
294 : std::vector<MathVector<WorldDim> > vGlobIP;
295 : std::vector<MathMatrix<dim, WorldDim> > vJT;
296 : std::vector<number> vValue;
297 :
298 : // iterate over all elements
299 0 : for(; iter != iterEnd; ++iter)
300 : {
301 : // get element
302 : grid_base_object* pElem = *iter;
303 :
304 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
305 0 : ReferenceObjectID roid = (ReferenceObjectID) pElem->reference_object_id();
306 :
307 : // get quadrature Rule for reference object id and order
308 : try{
309 : const QuadratureRule<dim>& rQuadRule
310 0 : = QuadratureRuleProvider<dim>::get(roid, quadOrder, type);
311 :
312 : // get reference element mapping by reference object id
313 : DimReferenceMapping<dim, WorldDim>& mapping
314 0 : = ReferenceMappingProvider::get<dim, WorldDim>(roid);
315 :
316 : // number of integration points
317 : const size_t numIP = rQuadRule.size();
318 :
319 : // get all corner coordinates
320 : CollectCornerCoordinates(vCorner, *pElem, aaPos, true);
321 :
322 : // update the reference mapping for the corners
323 0 : mapping.update(&vCorner[0]);
324 :
325 : // compute global integration points
326 0 : vGlobIP.resize(numIP);
327 0 : mapping.local_to_global(&(vGlobIP[0]), rQuadRule.points(), numIP);
328 :
329 : // compute transformation matrices
330 0 : vJT.resize(numIP);
331 0 : mapping.jacobian_transposed(&(vJT[0]), rQuadRule.points(), numIP);
332 :
333 : // compute integrand values at integration points
334 0 : vValue.resize(numIP);
335 : try
336 : {
337 0 : integrand.values(&(vValue[0]), &(vGlobIP[0]),
338 : pElem, &vCorner[0], rQuadRule.points(),
339 : &(vJT[0]),
340 : numIP);
341 : }
342 0 : UG_CATCH_THROW("Unable to compute values of integrand at integration point.");
343 :
344 : // reset contribution of this element
345 : number intValElem = 0;
346 :
347 : // loop integration points
348 0 : for(size_t ip = 0; ip < numIP; ++ip)
349 : {
350 : // get quadrature weight
351 : const number weightIP = rQuadRule.weight(ip);
352 :
353 : // get determinate of mapping
354 0 : const number det = SqrtGramDeterminant(vJT[ip]);
355 :
356 : // add contribution of integration point
357 0 : intValElem += vValue[ip] * weightIP * det;
358 : }
359 :
360 : // add to global sum
361 0 : integral += intValElem;
362 0 : if(aaElemContribs.valid())
363 0 : aaElemContribs[pElem] = intValElem;
364 :
365 0 : }UG_CATCH_THROW("SumValuesOnElems failed.");
366 : } // end elem
367 :
368 : // return the summed integral contributions of all elements
369 0 : return integral;
370 0 : }
371 :
372 : template <typename TGridFunction, int dim>
373 0 : number IntegrateSubset(IIntegrand<number, TGridFunction::dim> &spIntegrand,
374 : TGridFunction& spGridFct,
375 : int si, int quadOrder, std::string quadType)
376 : {
377 : // integrate elements of subset
378 : typedef typename TGridFunction::template dim_traits<dim>::grid_base_object grid_base_object;
379 : typedef typename TGridFunction::template dim_traits<dim>::const_iterator const_iterator;
380 :
381 0 : spIntegrand.set_subset(si);
382 :
383 : return Integrate<TGridFunction::dim,dim,const_iterator>
384 0 : (spGridFct.template begin<grid_base_object>(si),
385 : spGridFct.template end<grid_base_object>(si),
386 0 : spGridFct.domain()->position_accessor(),
387 : spIntegrand,
388 0 : quadOrder, quadType);
389 : }
390 :
391 :
392 : template <typename TGridFunction>
393 0 : number IntegrateSubsets(IIntegrand<number, TGridFunction::dim> &spIntegrand,
394 : TGridFunction& spGridFct,
395 : const char* subsets, int quadOrder,
396 : std::string quadType = std::string())
397 : {
398 : // world dimensions
399 : static const int dim = TGridFunction::dim;
400 :
401 : // read subsets
402 0 : SubsetGroup ssGrp(spGridFct.domain()->subset_handler());
403 0 : if(subsets != NULL)
404 : {
405 0 : ssGrp.add(TokenizeString(subsets));
406 0 : UG_COND_THROW(!SameDimensionsInAllSubsets(ssGrp), "IntegrateSubsets: Subsets '"<<subsets<<"' do not have same dimension."
407 : "Cannot integrate on subsets of different dimensions.");
408 0 : UG_LOG("IntegrateSubsets for subsets="<<subsets<<"\n");
409 : }
410 : else
411 : {
412 : // add all subsets and remove lower dim subsets afterwards
413 0 : ssGrp.add_all();
414 0 : RemoveLowerDimSubsets(ssGrp);
415 : }
416 :
417 : // reset value
418 : number value = 0;
419 :
420 : // loop subsets
421 0 : for(size_t i = 0; i < ssGrp.size(); ++i)
422 : {
423 : // get subset index
424 : const int si = ssGrp[i];
425 :
426 : // check dimension
427 0 : UG_COND_THROW(ssGrp.dim(i) > dim, "IntegrateSubsets: Dimension of subset is "<<ssGrp.dim(i)<<", but "
428 : " world dimension is "<<dim<<". Cannot integrate this.");
429 :
430 : // integrate elements of subset
431 : try{
432 0 : switch(ssGrp.dim(i))
433 : {
434 : case DIM_SUBSET_EMPTY_GRID: break;
435 0 : case 1: value += IntegrateSubset<TGridFunction, 1>(spIntegrand, spGridFct, si, quadOrder, quadType); break;
436 0 : case 2: value += IntegrateSubset<TGridFunction, 2>(spIntegrand, spGridFct, si, quadOrder, quadType); break;
437 0 : case 3: value += IntegrateSubset<TGridFunction, 3>(spIntegrand, spGridFct, si, quadOrder, quadType); break;
438 0 : default: UG_THROW("IntegrateSubsets: Dimension "<<ssGrp.dim(i)<<" not supported. "
439 : " World dimension is "<<dim<<".");
440 : }
441 : }
442 0 : UG_CATCH_THROW("IntegrateSubsets: Integration failed on subset "<<si);
443 : }
444 :
445 : #ifdef UG_PARALLEL
446 : // sum over processes
447 : if(pcl::NumProcs() > 1)
448 : {
449 : pcl::ProcessCommunicator com;
450 : number local = value;
451 : com.allreduce(&local, &value, 1, PCL_DT_DOUBLE, PCL_RO_SUM);
452 : }
453 : #endif
454 :
455 : // return the result
456 0 : return value;
457 : }
458 :
459 :
460 : ////////////////////////////////////////////////////////////////////////////////
461 : // UserData Integrand
462 : ////////////////////////////////////////////////////////////////////////////////
463 :
464 : //! For arbitrary UserData \f$\rho\f$, this class defines the integrand \f$\rho(u)\f$.
465 : template <typename TData, typename TGridFunction>
466 0 : class UserDataIntegrand
467 : : public StdIntegrand<TData, TGridFunction::dim, UserDataIntegrand<TData, TGridFunction> >
468 : {
469 : public:
470 : // world dimension of grid function
471 : static const int worldDim = TGridFunction::dim;
472 :
473 : // data type
474 : typedef TData data_type;
475 :
476 : private:
477 : // data to integrate
478 : SmartPtr<UserData<TData, worldDim> > m_spData;
479 :
480 : // grid function
481 : TGridFunction* m_spGridFct;
482 :
483 : // time
484 : number m_time;
485 :
486 : public:
487 : /// constructor
488 0 : UserDataIntegrand(SmartPtr<UserData<TData, worldDim> > spData,
489 : TGridFunction* spGridFct,
490 : number time)
491 0 : : m_spData(spData), m_spGridFct(spGridFct), m_time(time)
492 : {
493 0 : m_spData->set_function_pattern(spGridFct->function_pattern());
494 0 : };
495 :
496 : /// constructor
497 : UserDataIntegrand(SmartPtr<UserData<TData, worldDim> > spData,
498 : number time)
499 : : m_spData(spData), m_spGridFct(NULL), m_time(time)
500 : {
501 : UG_COND_THROW(m_spData->requires_grid_fct(),
502 : "UserDataIntegrand: Missing GridFunction, but data requires grid function.");
503 : };
504 :
505 : /// \copydoc IIntegrand::values
506 : template <int elemDim>
507 0 : void evaluate(TData vValue[],
508 : const MathVector<worldDim> vGlobIP[],
509 : GridObject* pElem,
510 : const MathVector<worldDim> vCornerCoords[],
511 : const MathVector<elemDim> vLocIP[],
512 : const MathMatrix<elemDim, worldDim> vJT[],
513 : const size_t numIP)
514 : {
515 : // get local solution if needed
516 0 : if(m_spData->requires_grid_fct())
517 : {
518 : // create storage
519 : LocalIndices ind;
520 : LocalVector u;
521 :
522 : // get global indices
523 0 : m_spGridFct->indices(pElem, ind);
524 :
525 : // adapt local algebra
526 0 : u.resize(ind);
527 :
528 : // read local values of u
529 0 : GetLocalVector(u, *m_spGridFct);
530 :
531 : // compute data
532 : try{
533 0 : (*m_spData)(vValue, vGlobIP, m_time, this->m_si, pElem,
534 : vCornerCoords, vLocIP, numIP, &u, &vJT[0]);
535 : }
536 0 : UG_CATCH_THROW("UserDataIntegrand: Cannot evaluate data.");
537 : }
538 : else
539 : {
540 : // compute data
541 : try{
542 0 : (*m_spData)(vValue, vGlobIP, m_time, this->m_si, numIP);
543 : }
544 0 : UG_CATCH_THROW("UserDataIntegrand: Cannot evaluate data.");
545 : }
546 :
547 0 : }
548 : };
549 :
550 :
551 :
552 :
553 : //! For arbitrary UserData \f$f\f$ (of type TData), this class defines the integrand \f$f^2(u)\f$.
554 : template <typename TData, typename TGridFunction>
555 0 : class UserDataIntegrandSq
556 : : public StdIntegrand<number, TGridFunction::dim, UserDataIntegrandSq<TData, TGridFunction> >
557 : {
558 : public:
559 : // world dimension of grid function
560 : static const int worldDim = TGridFunction::dim;
561 :
562 : // data type
563 : typedef TData data_type;
564 :
565 : private:
566 : // data to integrate
567 : SmartPtr<UserData<TData, worldDim> > m_spData;
568 :
569 : // grid function
570 : const TGridFunction* m_pGridFct;
571 :
572 : // time
573 : number m_time;
574 :
575 : public:
576 : /// constructor
577 0 : UserDataIntegrandSq(SmartPtr<UserData<TData, worldDim> > spData,
578 : const TGridFunction* pGridFct,
579 : number time)
580 0 : : m_spData(spData), m_pGridFct(pGridFct), m_time(time)
581 : {
582 0 : m_spData->set_function_pattern(pGridFct->function_pattern());
583 0 : };
584 :
585 : /// constructor
586 : UserDataIntegrandSq(SmartPtr<UserData<TData, worldDim> > spData,
587 : number time)
588 : : m_spData(spData), m_pGridFct(NULL), m_time(time)
589 : {
590 : if(m_spData->requires_grid_fct())
591 : UG_THROW("UserDataIntegrand: Missing GridFunction, but "
592 : " data requires grid function.")
593 : };
594 :
595 :
596 : /// \copydoc IIntegrand::values
597 : template <int elemDim>
598 0 : void evaluate(number vValue[],
599 : const MathVector<worldDim> vGlobIP[],
600 : GridObject* pElem,
601 : const MathVector<worldDim> vCornerCoords[],
602 : const MathVector<elemDim> vLocIP[],
603 : const MathMatrix<elemDim, worldDim> vJT[],
604 : const size_t numIP)
605 : {
606 :
607 0 : std::vector<TData> tmpValues(numIP);
608 :
609 : // get local solution if needed
610 0 : if(m_spData->requires_grid_fct())
611 : {
612 : UG_ASSERT(m_pGridFct!=NULL, "Error: Requires valid pointer!")
613 : // create storage
614 : LocalIndices ind;
615 : LocalVector u;
616 :
617 : // get global indices
618 0 : m_pGridFct->indices(pElem, ind);
619 :
620 : // adapt local algebra
621 0 : u.resize(ind);
622 :
623 : // read local values of u
624 0 : GetLocalVector(u, *m_pGridFct);
625 :
626 : // compute data
627 : try{
628 0 : (*m_spData)(&tmpValues.front(), vGlobIP, m_time, this->m_si, pElem,
629 : vCornerCoords, vLocIP, numIP, &u, &vJT[0]);
630 : }
631 0 : UG_CATCH_THROW("UserDataIntegrand: Cannot evaluate data.");
632 : }
633 : else
634 : {
635 : // compute data
636 : try{
637 0 : (*m_spData)(&tmpValues.front(), vGlobIP, m_time, this->m_si, numIP);
638 : }
639 0 : UG_CATCH_THROW("UserDataIntegrand: Cannot evaluate data.");
640 : }
641 :
642 0 : for (size_t i=0; i<numIP; ++i)
643 : {
644 0 : vValue[i]=inner_prod(tmpValues[i], tmpValues[i]);
645 : }
646 :
647 0 : }
648 : protected:
649 :
650 : number inner_prod(const number &d1, const number &d2)
651 : {return d1*d2;}
652 :
653 : number inner_prod(const MathVector<worldDim> &d1, const MathVector<worldDim> &d2)
654 : {return VecDot(d1, d2);}
655 :
656 : template<typename T>
657 : number inner_prod(const T &d1, const T &d2)
658 : { UG_ASSERT(0, "NOT IMPLEMENTED"); return 0.0;}
659 : };
660 :
661 :
662 :
663 :
664 : //! For arbitrary UserData \f$f\f$ and grid functions \f$u_1\f$ and \f$u_2\f$, this class (should) define the integrand \f$ (f(u_1)- f(u_2))^2 \f$.
665 : template <typename TData, typename TGridFunction>
666 : class UserDataDistIntegrandSq
667 : : public StdIntegrand<number, TGridFunction::dim, UserDataDistIntegrandSq<TData, TGridFunction> >
668 : {
669 : public:
670 : /// world dimension of grid function
671 : static const int worldDim = TGridFunction::dim;
672 : // typedef typename L2Integrand<TGridFunction>::weight_type weight_type;
673 :
674 : protected:
675 : ScalarGridFunctionData<TGridFunction> m_fineData;
676 : const int m_fineTopLevel;
677 :
678 : ScalarGridFunctionData<TGridFunction> m_coarseData;
679 : const int m_coarseTopLevel;
680 :
681 : /// multigrid
682 : SmartPtr<MultiGrid> m_spMG;
683 :
684 :
685 :
686 : SmartPtr<UserData<TData, worldDim> > m_spData;
687 : // ConstSmartPtr<weight_type> m_spWeight;
688 :
689 : double m_time;
690 : public:
691 :
692 : /// constructor (1st is fine grid function)
693 : UserDataDistIntegrandSq(SmartPtr<UserData<TData, worldDim> > spData, TGridFunction& fineGridFct, size_t fineCmp,
694 : TGridFunction& coarseGridFct, size_t coarseCmp)
695 : : m_fineData(fineGridFct, fineCmp), m_fineTopLevel(fineGridFct.dof_distribution()->grid_level().level()),
696 : m_coarseData(coarseGridFct, coarseCmp), m_coarseTopLevel(coarseGridFct.dof_distribution()->grid_level().level()),
697 : m_spMG(m_fineData.domain()->grid()), m_spData(spData), m_time(0.0) /*, m_spWeight(make_sp(new ConstUserNumber<TGridFunction::dim>(1.0)))*/
698 : {
699 : if(m_fineTopLevel < m_coarseTopLevel)
700 : UG_THROW("UserDataDistIntegrandSq: fine and top level inverted.");
701 :
702 : if(m_fineData.domain().get() != m_coarseData.domain().get())
703 : UG_THROW("UserDataDistIntegrandSq: grid functions defined on different domains.");
704 : };
705 :
706 : virtual ~UserDataDistIntegrandSq() {}
707 :
708 : /// sets subset
709 : virtual void set_subset(int si)
710 : {
711 : if(!m_fineData.is_def_in_subset(si))
712 : UG_THROW("UserDataDistIntegrandSq: Grid function component"
713 : <<m_fineData.fct()<<" not defined on subset "<<si);
714 : if(!m_coarseData.is_def_in_subset(si))
715 : UG_THROW("UserDataDistIntegrandSq: Grid function component"
716 : <<m_coarseData.fct()<<" not defined on subset "<<si);
717 : IIntegrand<number, worldDim>::set_subset(si);
718 : }
719 :
720 : /// \copydoc IIntegrand::values
721 : template <int elemDim>
722 : void evaluate(number vValue[],
723 : const MathVector<worldDim> vGlobIP[],
724 : GridObject* pFineElem,
725 : const MathVector<worldDim> vCornerCoords[],
726 : const MathVector<elemDim> vFineLocIP[],
727 : const MathMatrix<elemDim, worldDim> vJT[],
728 : const size_t numIP)
729 : {
730 :
731 :
732 : // must return 0.0, if m_spData is independent of grid functions
733 : if(!m_spData->requires_grid_fct()) {
734 : for (size_t i=0; i<numIP; ++i) { vValue[i]=0.0; }
735 : return;
736 : }
737 :
738 : // get coarse element
739 : GridObject* pCoarseElem = pFineElem;
740 : if(m_coarseTopLevel < m_fineTopLevel){
741 : int parentLevel = m_spMG->get_level(pCoarseElem);
742 : while(parentLevel > m_coarseTopLevel){
743 : pCoarseElem = m_spMG->get_parent(pCoarseElem);
744 : parentLevel = m_spMG->get_level(pCoarseElem);
745 : }
746 : }
747 :
748 : // get corner coordinates
749 : typedef typename TGridFunction::template dim_traits<elemDim>::grid_base_object TElem;
750 : std::vector<MathVector<worldDim> > vCornerCoarse;
751 : CollectCornerCoordinates(vCornerCoarse, *static_cast<TElem*>(pCoarseElem), *m_coarseData.domain());
752 :
753 : // get Reference Mapping
754 : const ReferenceObjectID coarseROID = pCoarseElem->reference_object_id();
755 : DimReferenceMapping<elemDim, worldDim>& map
756 : = ReferenceMappingProvider::get<elemDim, worldDim>(coarseROID, vCornerCoarse);
757 :
758 : std::vector<MathVector<elemDim> > vCoarseLocIP;
759 : vCoarseLocIP.resize(numIP);
760 : for(size_t ip = 0; ip < vCoarseLocIP.size(); ++ip) VecSet(vCoarseLocIP[ip], 0.0);
761 : map.global_to_local(&vCoarseLocIP[0], vGlobIP, numIP);
762 :
763 : // element weights
764 : /* typedef typename weight_type::data_type ipdata_type;
765 : std::vector<ipdata_type> fineElemWeights(numIP, 1.0);
766 : UG_ASSERT(m_spWeight.valid(), "L2DistIntegrand::evaluate requires valid weights! ");
767 : (*m_spWeight)(&fineElemWeights[0], vFineGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), numIP);*/
768 :
769 :
770 : // get trial space
771 : /* const LocalShapeFunctionSet<elemDim>& rFineLSFS =
772 : LocalFiniteElementProvider::get<elemDim>(fineROID, m_fineData.id());
773 : const LocalShapeFunctionSet<elemDim>& rCoarseLSFS =
774 : LocalFiniteElementProvider::get<elemDim>(coarseROID, m_coarseData.id());*/
775 :
776 : // get multiindices of element
777 : /*std::vector<DoFIndex> vFineMI, vCoarseMI;
778 : m_fineData.dof_indices(pFineElem, vFineMI);
779 : m_coarseData.dof_indices(pCoarseElem, vCoarseMI);*/
780 :
781 : TData fineValues[numIP];
782 : TData coarseValues[numIP];
783 :
784 :
785 :
786 : // compute coarse data
787 : try{
788 : LocalVector uCoarse;
789 : LocalIndices indCoarse;
790 :
791 : m_coarseData.grid_function().indices(pCoarseElem, indCoarse);
792 : uCoarse.resize(indCoarse);
793 : GetLocalVector(uCoarse, m_coarseData.grid_function());
794 :
795 : (*m_spData)(coarseValues, vGlobIP, m_time, this->m_si, pCoarseElem,
796 : &vCornerCoarse[0], &vCoarseLocIP[0], numIP, &uCoarse, &vJT[0]);
797 : } UG_CATCH_THROW("UserDataDistIntegrandSq: Cannot evaluate coarse data.");
798 :
799 :
800 : // compute fine data
801 : try{
802 : LocalVector uFine;
803 : LocalIndices indFine;
804 :
805 : m_fineData.grid_function().indices(pFineElem, indFine);
806 : uFine.resize(indFine);
807 : GetLocalVector(uFine, m_fineData.grid_function());
808 :
809 : (*m_spData)(fineValues, vGlobIP, m_time, this->m_si, pFineElem,
810 : vCornerCoords, vFineLocIP, numIP, &uFine, &vJT[0]);
811 : } UG_CATCH_THROW("UserDataDistIntegrandSq: Cannot evaluate fine data.");
812 :
813 : // loop all integration points
814 : for(size_t ip = 0; ip < numIP; ++ip)
815 : {
816 : vValue[ip] = inner_dist2(fineValues[ip], coarseValues[ip]);
817 : }
818 :
819 : }
820 :
821 :
822 : protected:
823 :
824 : number inner_prod(const number &d1, const number &d2)
825 : {return d1*d2;}
826 :
827 : number inner_prod(const MathVector<worldDim> &d1, const MathVector<worldDim> &d2)
828 : {return VecDot(d1, d2);}
829 :
830 : template<typename T>
831 : number inner_prod(const T &d1, const T &d2)
832 : { UG_ASSERT(0, "NOT IMPLEMENTED"); return 0.0;}
833 :
834 :
835 :
836 : number inner_dist2(const number &v1, const number &v2)
837 : { return (v1-v2)*(v1-v2); }
838 :
839 : number inner_dist2(const MathVector<worldDim> &v1, const MathVector<worldDim> &v2)
840 : { return VecDistanceSq(v1, v2); }
841 :
842 : template<typename T>
843 : number inner_dist2(const T &d1, const T &d2)
844 : { UG_ASSERT(0, "NOT IMPLEMENTED"); return 0.0;}
845 :
846 : };
847 :
848 :
849 : ////////////////////////////////////////////////////////////////////////////////
850 : ////////////////////////////////////////////////////////////////////////////////
851 : // Volume Integrand implementations
852 : ////////////////////////////////////////////////////////////////////////////////
853 : ////////////////////////////////////////////////////////////////////////////////
854 :
855 : template <typename TGridFunction>
856 0 : number Integral(SmartPtr<UserData<number, TGridFunction::dim> > spData,
857 : TGridFunction& spGridFct,
858 : const char* subsets, number time,
859 : int quadOrder, std::string quadType)
860 : {
861 0 : UserDataIntegrand<number, TGridFunction> spIntegrand(spData, &spGridFct, time);
862 0 : return IntegrateSubsets(spIntegrand, spGridFct, subsets, quadOrder, quadType);
863 : }
864 :
865 : template <typename TGridFunction>
866 0 : number Integral(SmartPtr<UserData<number, TGridFunction::dim> > spData, SmartPtr<TGridFunction> spGridFct,
867 : const char* subsets, number time, int quadOrder, std::string quadType)
868 0 : { return Integral(spData, *spGridFct, subsets, time, quadOrder, quadType); }
869 :
870 : template <typename TGridFunction>
871 0 : number Integral(SmartPtr<UserData<number, TGridFunction::dim> > spData, SmartPtr<TGridFunction> spGridFct,const char* subsets,number time, int order)
872 0 : {return Integral(spData, spGridFct, subsets, time, order, "best");}
873 :
874 : template <typename TGridFunction>
875 0 : number Integral(SmartPtr<UserData<number, TGridFunction::dim> > spData, SmartPtr<TGridFunction> spGridFct,const char* subsets,number time)
876 0 : {return Integral(spData, spGridFct, subsets, time, 1, "best");}
877 :
878 : template <typename TGridFunction>
879 0 : number Integral(SmartPtr<UserData<number, TGridFunction::dim> > spData, SmartPtr<TGridFunction> spGridFct,number time)
880 0 : {return Integral(spData, spGridFct, NULL, time, 1, "best");}
881 :
882 : template <typename TGridFunction>
883 0 : number Integral(SmartPtr<UserData<number, TGridFunction::dim> > spData, SmartPtr<TGridFunction> spGridFct,const char* subsets)
884 0 : {return Integral(spData, spGridFct, subsets, 0.0, 1, "best");}
885 :
886 : template <typename TGridFunction>
887 0 : number Integral(SmartPtr<UserData<number, TGridFunction::dim> > spData, SmartPtr<TGridFunction> spGridFct)
888 0 : {return Integral(spData, spGridFct, NULL, 0.0, 1, "best");}
889 :
890 : ///////////////
891 : // const data
892 : ///////////////
893 :
894 : template <typename TGridFunction>
895 0 : number Integral(number val, SmartPtr<TGridFunction> spGridFct,
896 : const char* subsets,
897 : number time, int quadOrder)
898 : {
899 : static const int dim = TGridFunction::dim;
900 0 : SmartPtr<UserData<number, dim> > sp =
901 0 : make_sp(new ConstUserNumber<dim>(val));
902 0 : return Integral(sp, spGridFct, subsets, time, quadOrder);
903 : }
904 :
905 : template <typename TGridFunction>
906 0 : number Integral(number val, SmartPtr<TGridFunction> spGridFct,const char* subsets,number time)
907 0 : {return Integral(val, spGridFct, subsets, time, 1);}
908 :
909 : template <typename TGridFunction>
910 0 : number Integral(number val, SmartPtr<TGridFunction> spGridFct,number time)
911 0 : {return Integral(val, spGridFct, NULL, time, 1);}
912 :
913 : template <typename TGridFunction>
914 0 : number Integral(number val, SmartPtr<TGridFunction> spGridFct,const char* subsets)
915 0 : {return Integral(val, spGridFct, subsets, 0.0, 1);}
916 :
917 : template <typename TGridFunction>
918 0 : number Integral(number val, SmartPtr<TGridFunction> spGridFct)
919 0 : {return Integral(val, spGridFct, NULL, 0.0, 1);}
920 :
921 : ///////////////
922 : // lua data
923 : ///////////////
924 :
925 : #ifdef UG_FOR_LUA
926 : template <typename TGridFunction>
927 0 : number Integral(const char* luaFct,
928 : SmartPtr<TGridFunction> spGridFct,
929 : const char* subsets, number time,
930 : int quadOrder, std::string quadType)
931 : {
932 : static const int dim = TGridFunction::dim;
933 0 : SmartPtr<UserData<number, dim> > sp =
934 : LuaUserDataFactory<number, dim>::create(luaFct);
935 0 : return Integral(sp, *spGridFct, subsets, time, quadOrder, quadType);
936 : }
937 :
938 : template <typename TGridFunction>
939 0 : number Integral(const char* luaFct, SmartPtr<TGridFunction> spGridFct,const char* subsets,number time, int quadOrder)
940 0 : {return Integral(luaFct, spGridFct, subsets, time, quadOrder, "best");}
941 :
942 : template <typename TGridFunction>
943 0 : number Integral(const char* luaFct, SmartPtr<TGridFunction> spGridFct,const char* subsets,number time)
944 0 : {return Integral(luaFct, spGridFct, subsets, time, 1, "best");}
945 :
946 : template <typename TGridFunction>
947 0 : number Integral(const char* luaFct, SmartPtr<TGridFunction> spGridFct,number time)
948 0 : {return Integral(luaFct, spGridFct, NULL, time, 1, "best");}
949 :
950 : template <typename TGridFunction>
951 0 : number Integral(const char* luaFct, SmartPtr<TGridFunction> spGridFct,const char* subsets)
952 0 : {return Integral(luaFct, spGridFct, subsets, 0.0, 1, "best");}
953 :
954 : template <typename TGridFunction>
955 0 : number Integral(const char* luaFct, SmartPtr<TGridFunction> spGridFct)
956 0 : {return Integral(luaFct, spGridFct, NULL, 0.0, 1, "best");}
957 :
958 : #endif
959 :
960 :
961 :
962 : template <typename TGridFunction>
963 : class MaximumDistIntegrand
964 : : public StdIntegrand<number, TGridFunction::dim, MaximumDistIntegrand<TGridFunction> >
965 : {
966 : public:
967 : /// world dimension of grid function
968 : static const int worldDim = TGridFunction::dim;
969 :
970 : protected:
971 : ScalarGridFunctionData<TGridFunction> m_scalarData;
972 :
973 : double m_max, m_min;
974 :
975 : public:
976 : /// constructor
977 0 : MaximumDistIntegrand(TGridFunction& gridFct, size_t cmp)
978 :
979 : : m_scalarData(gridFct, cmp),
980 0 : m_max(-std::numeric_limits<float>::infinity()),
981 0 : m_min(std::numeric_limits<float>::infinity())
982 : {};
983 :
984 0 : virtual ~MaximumDistIntegrand() {};
985 :
986 : /// sets subset
987 0 : virtual void set_subset(int si)
988 : {
989 0 : UG_COND_THROW(!m_scalarData.is_def_in_subset(si), "L2ErrorIntegrand: Grid function component" <<
990 : m_scalarData.fct()<<" not defined on subset "<<si);
991 : IIntegrand<number, worldDim>::set_subset(si);
992 0 : }
993 :
994 : void reset()
995 : {
996 : m_min = std::numeric_limits<double>::infinity();
997 : m_max = -std::numeric_limits<double>::infinity();
998 : }
999 :
1000 0 : double min() { return m_min; }
1001 : double max() { return m_max; }
1002 :
1003 : /// \copydoc IIntegrand::values
1004 : template <int elemDim>
1005 0 : void evaluate(number vValue[],
1006 : const MathVector<worldDim> vGlobIP[],
1007 : GridObject* pElem,
1008 : const MathVector<worldDim> vCornerCoords[],
1009 : const MathVector<elemDim> vLocIP[],
1010 : const MathMatrix<elemDim, worldDim> vJT[],
1011 : const size_t numIP)
1012 : {
1013 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
1014 0 : const ReferenceObjectID roid = pElem->reference_object_id();
1015 :
1016 : try{
1017 : // get trial space
1018 : const LocalShapeFunctionSet<elemDim>& rTrialSpace =
1019 : LocalFiniteElementProvider::get<elemDim>(roid, m_scalarData.id());
1020 :
1021 : // number of dofs on element
1022 0 : const size_t num_sh = rTrialSpace.num_sh();
1023 :
1024 : // get multiindices of element
1025 : std::vector<DoFIndex> ind; // aux. index array
1026 : m_scalarData.dof_indices(pElem, ind);
1027 :
1028 : // check multi indices
1029 0 : UG_COND_THROW(ind.size() != num_sh, "L2ErrorIntegrand::evaluate: Wrong number of multi indices.");
1030 :
1031 : // evaluate function in corners.
1032 0 : for(size_t sh = 0; sh < num_sh; ++sh)
1033 : {
1034 : // get value at shape point (e.g. corner for P1 fct)
1035 0 : const number val = DoFRef(m_scalarData.grid_function(), ind[sh]);
1036 :
1037 0 : m_max = std::max(m_max, val);
1038 0 : m_min = std::min(m_min, val);
1039 : }
1040 :
1041 :
1042 0 : } UG_CATCH_THROW("L2ErrorIntegrand::evaluate: trial space missing.");
1043 0 : }
1044 : };
1045 :
1046 :
1047 : template <typename TGridFunction>
1048 0 : number GetMinimum(TGridFunction &gridFct, const char* cmp, const char* subsets)
1049 : {
1050 : // get function id of name & check that function exists
1051 : const size_t fct = gridFct.fct_id_by_name(cmp);
1052 0 : UG_COND_THROW(fct >= gridFct.num_fct(),
1053 : "L2Error: Function space does not contain a function with name " << cmp << ".");
1054 :
1055 : // Find minimum.
1056 : MaximumDistIntegrand<TGridFunction> spIntegrand(gridFct, fct);
1057 0 : IntegrateSubsets(spIntegrand, gridFct, subsets, 2);
1058 0 : return spIntegrand.min();
1059 : }
1060 :
1061 : template <typename TGridFunction>
1062 0 : number Minimum(SmartPtr<TGridFunction> spGridFct, const char* cmp, const char* subsets)
1063 0 : { return GetMinimum(*spGridFct, cmp, subsets); }
1064 :
1065 :
1066 : ////////////////////////////////////////////////////////////////////////////////
1067 : // L2 Error Integrand
1068 : ////////////////////////////////////////////////////////////////////////////////
1069 :
1070 : template <typename TGridFunction>
1071 : class L2ErrorIntegrand
1072 : : public StdIntegrand<number, TGridFunction::dim, L2ErrorIntegrand<TGridFunction> >
1073 : {
1074 : public:
1075 : /// world dimension of grid function
1076 : static const int worldDim = TGridFunction::dim;
1077 :
1078 : private:
1079 : ScalarGridFunctionData<TGridFunction> m_scalarData;
1080 :
1081 : /// exact solution
1082 : SmartPtr<UserData<number, worldDim> > m_spExactSolution;
1083 :
1084 : /// time
1085 : number m_time;
1086 :
1087 : public:
1088 : /// constructor
1089 0 : L2ErrorIntegrand(SmartPtr<UserData<number, worldDim> > spExactSol,
1090 : TGridFunction& gridFct, size_t cmp,
1091 : number time)
1092 : : m_scalarData(gridFct, cmp),
1093 0 : m_spExactSolution(spExactSol), m_time(time)
1094 : {};
1095 :
1096 0 : virtual ~L2ErrorIntegrand() {};
1097 :
1098 : /// sets subset
1099 0 : virtual void set_subset(int si)
1100 : {
1101 0 : if(!m_scalarData.is_def_in_subset(si))
1102 0 : UG_THROW("L2ErrorIntegrand: Grid function component"
1103 : <<m_scalarData.fct()<<" not defined on subset "<<si);
1104 : IIntegrand<number, worldDim>::set_subset(si);
1105 0 : }
1106 :
1107 : /// \copydoc IIntegrand::values
1108 : template <int elemDim>
1109 0 : void evaluate(number vValue[],
1110 : const MathVector<worldDim> vGlobIP[],
1111 : GridObject* pElem,
1112 : const MathVector<worldDim> vCornerCoords[],
1113 : const MathVector<elemDim> vLocIP[],
1114 : const MathMatrix<elemDim, worldDim> vJT[],
1115 : const size_t numIP)
1116 : {
1117 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
1118 0 : const ReferenceObjectID roid = pElem->reference_object_id();
1119 :
1120 : try{
1121 : // get trial space
1122 : const LocalShapeFunctionSet<elemDim>& rTrialSpace =
1123 : LocalFiniteElementProvider::get<elemDim>(roid, m_scalarData.id());
1124 :
1125 : // number of dofs on element
1126 0 : const size_t num_sh = rTrialSpace.num_sh();
1127 :
1128 : // get multiindices of element
1129 : std::vector<DoFIndex> ind; // aux. index array
1130 : m_scalarData.dof_indices(pElem, ind);
1131 :
1132 : // check multi indices
1133 0 : if(ind.size() != num_sh)
1134 0 : UG_THROW("L2ErrorIntegrand::evaluate: Wrong number of"
1135 : " multi indices.");
1136 :
1137 : // loop all integration points
1138 0 : for(size_t ip = 0; ip < numIP; ++ip)
1139 : {
1140 : // compute exact solution at integration point
1141 : number exactSolIP;
1142 0 : (*m_spExactSolution)(exactSolIP, vGlobIP[ip], m_time, this->subset());
1143 :
1144 : // compute approximated solution at integration point
1145 : number approxSolIP = 0.0;
1146 0 : for(size_t sh = 0; sh < num_sh; ++sh)
1147 : {
1148 : // get value at shape point (e.g. corner for P1 fct)
1149 0 : const number valSH = DoFRef(m_scalarData.grid_function(), ind[sh]);
1150 :
1151 : // add shape fct at ip * value at shape
1152 0 : approxSolIP += valSH * rTrialSpace.shape(sh, vLocIP[ip]);
1153 : }
1154 :
1155 : // get squared of difference
1156 0 : vValue[ip] = (exactSolIP - approxSolIP);
1157 0 : vValue[ip] *= vValue[ip];
1158 : }
1159 :
1160 : }
1161 0 : UG_CATCH_THROW("L2ErrorIntegrand::evaluate: trial space missing.");
1162 0 : }
1163 : };
1164 :
1165 :
1166 :
1167 : ////////////////////////////////////////////////////////////////////////////////
1168 : // L2 Error Integrand
1169 : ////////////////////////////////////////////////////////////////////////////////
1170 :
1171 : /// computes the l2 error function on the whole domain or on some subsets
1172 : /**
1173 : * This function computes the L2-difference between a given analytic function
1174 : * and a grid function.
1175 : *
1176 : * \param[in] ExactSol analytic function
1177 : * \param[in] spGridFct grid function
1178 : * \param[in] cmp symbolic name of component function
1179 : * \param[in] time time point
1180 : * \param[in] quadOrder order of quadrature rule
1181 : * \param[in] subsets subsets, where to compute
1182 : * \returns number l2-norm of difference
1183 : */
1184 : template <typename TGridFunction>
1185 0 : number L2Error(SmartPtr<UserData<number, TGridFunction::dim> > spExactSol,
1186 : TGridFunction& gridFct, const char* cmp,
1187 : number time, int quadOrder, const char* subsets)
1188 : {
1189 : // get function id of name
1190 : const size_t fct = gridFct.fct_id_by_name(cmp);
1191 :
1192 : // check that function exists
1193 0 : UG_COND_THROW(fct >= gridFct.num_fct(),
1194 : "L2Error: Function space does not contain a function with name " << cmp << ".");
1195 :
1196 0 : L2ErrorIntegrand<TGridFunction> spIntegrand(spExactSol, gridFct, fct, time);
1197 0 : return sqrt(IntegrateSubsets(spIntegrand, gridFct, subsets, quadOrder));
1198 : }
1199 :
1200 : template <typename TGridFunction>
1201 0 : number L2Error(SmartPtr<UserData<number, TGridFunction::dim> > spExactSol,
1202 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
1203 : number time, int quadOrder, const char* subsets)
1204 0 : { return L2Error(spExactSol, *spGridFct, cmp, time, quadOrder, subsets); }
1205 :
1206 : template <typename TGridFunction>
1207 0 : number L2Error(SmartPtr<UserData<number, TGridFunction::dim> > spExactSol,
1208 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
1209 : number time, int quadOrder)
1210 0 : { return L2Error(spExactSol, *spGridFct, cmp, time, quadOrder, NULL); }
1211 :
1212 :
1213 :
1214 :
1215 : #ifdef UG_FOR_LUA
1216 : template <typename TGridFunction>
1217 0 : number L2Error(const char* ExactSol,
1218 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
1219 : number time, int quadOrder, const char* subsets)
1220 : {
1221 0 : SmartPtr<UserData<number, TGridFunction::dim> > spExactSol
1222 0 : = make_sp(new LuaUserData<number, TGridFunction::domain_type::dim>(ExactSol));
1223 0 : return L2Error(spExactSol, spGridFct, cmp, time, quadOrder, subsets);
1224 : }
1225 :
1226 : template <typename TGridFunction>
1227 0 : number L2Error(const char* ExactSol,
1228 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
1229 : number time, int quadOrder)
1230 : {
1231 0 : return L2Error(ExactSol, spGridFct, cmp, time, quadOrder, NULL);
1232 : }
1233 : #endif //UG_FOR_LUA
1234 :
1235 :
1236 :
1237 : #ifdef UG_DEBUG
1238 : // True, if max norm is less than SMALL.
1239 : template <int worldDim, int elemDim>
1240 : void CheckGeneralizedInverse(const MathMatrix<elemDim, worldDim> &JT, const MathMatrix<worldDim, elemDim> &JTInv)
1241 : {
1242 : if (worldDim==elemDim) return;
1243 :
1244 : MathMatrix<worldDim, elemDim> myTmp;
1245 : GeneralizedInverse(myTmp, JT); //old algorithms with Inverse, but Inverse is only defined for NxN Matrix
1246 : MatSubtract(myTmp, JTInv, myTmp);
1247 :
1248 : UG_ASSERT(MatMaxNorm(myTmp)<SMALL,
1249 : "The inverse matrix is not identity to the map jacobian transposed inverse.");
1250 :
1251 : }
1252 : #endif
1253 :
1254 :
1255 : ////////////////////////////////////////////////////////////////////////////////
1256 : // H1 Error Integrand
1257 : ////////////////////////////////////////////////////////////////////////////////
1258 :
1259 : template <typename TGridFunction>
1260 0 : class H1ErrorIntegrand
1261 : : public StdIntegrand<number, TGridFunction::dim, H1ErrorIntegrand<TGridFunction> >
1262 : {
1263 : public:
1264 : /// world dimension of grid function
1265 : static const int worldDim = TGridFunction::dim;
1266 :
1267 : private:
1268 : /// grid function
1269 : ScalarGridFunctionData<TGridFunction> m_scalarData;
1270 :
1271 : /// exact solution
1272 : SmartPtr<UserData<number, worldDim> > m_spExactSolution;
1273 :
1274 : /// exact gradient
1275 : SmartPtr<UserData<MathVector<worldDim>, worldDim> > m_spExactGrad;
1276 :
1277 : /// time
1278 : number m_time;
1279 :
1280 : public:
1281 : /// constructor
1282 0 : H1ErrorIntegrand(SmartPtr<UserData<number, worldDim> > spExactSol,
1283 : SmartPtr<UserData<MathVector<worldDim>, worldDim> > spExactGrad,
1284 : TGridFunction& gridFct, size_t cmp,
1285 : number time)
1286 : : m_scalarData (gridFct, cmp),
1287 : m_spExactSolution(spExactSol),
1288 : m_spExactGrad(spExactGrad),
1289 0 : m_time(time)
1290 0 : {}
1291 :
1292 : /// sets subset
1293 0 : virtual void set_subset(int si)
1294 : {
1295 0 : if(!m_scalarData.is_def_in_subset(si))
1296 0 : UG_THROW("H1Error: Grid function component"
1297 : <<m_scalarData.fct()<<" not defined on subset "<<si);
1298 : IIntegrand<number, worldDim>::set_subset(si);
1299 0 : }
1300 :
1301 : /// \copydoc IIntegrand::values
1302 : template <int elemDim>
1303 0 : void evaluate(number vValue[],
1304 : const MathVector<worldDim> vGlobIP[],
1305 : GridObject* pElem,
1306 : const MathVector<worldDim> vCornerCoords[],
1307 : const MathVector<elemDim> vLocIP[],
1308 : const MathMatrix<elemDim, worldDim> vJT[],
1309 : const size_t numIP)
1310 : {
1311 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
1312 0 : const ReferenceObjectID roid = pElem->reference_object_id();
1313 :
1314 : DimReferenceMapping<elemDim, worldDim>& map
1315 : = ReferenceMappingProvider::get<elemDim, worldDim>(roid, vCornerCoords);
1316 :
1317 : try{
1318 : // get trial space
1319 : const LocalShapeFunctionSet<elemDim>& rTrialSpace =
1320 : LocalFiniteElementProvider::get<elemDim>(roid, m_scalarData.id());
1321 :
1322 : // number of dofs on element
1323 0 : const size_t num_sh = rTrialSpace.num_sh();
1324 :
1325 : // get multiindices of element
1326 : std::vector<DoFIndex> ind; // aux. index array
1327 : m_scalarData.dof_indices(pElem, ind);
1328 :
1329 : // check multi indices
1330 0 : if(ind.size() != num_sh)
1331 0 : UG_THROW("H1ErrorIntegrand::evaluate: Wrong number of"
1332 : " multi indices.");
1333 :
1334 : // loop all integration points
1335 0 : std::vector<MathVector<elemDim> > vLocGradient(num_sh);
1336 0 : for(size_t ip = 0; ip < numIP; ++ip)
1337 : {
1338 : // compute exact solution at integration point
1339 : number exactSolIP;
1340 0 : (*m_spExactSolution)(exactSolIP, vGlobIP[ip], m_time, this->subset());
1341 :
1342 : // compute exact gradient at integration point
1343 : MathVector<worldDim> exactGradIP;
1344 0 : (*m_spExactGrad)(exactGradIP, vGlobIP[ip], m_time, this->subset());
1345 :
1346 : // compute shape gradients at ip
1347 0 : rTrialSpace.grads(&vLocGradient[0], vLocIP[ip]);
1348 :
1349 : // compute approximated solution at integration point
1350 : number approxSolIP = 0.0;
1351 : MathVector<elemDim> locTmp; VecSet(locTmp, 0.0);
1352 0 : for(size_t sh = 0; sh < num_sh; ++sh)
1353 : {
1354 : // get value at shape point (e.g. corner for P1 fct)
1355 0 : const number valSH = DoFRef(m_scalarData.grid_function(), ind[sh]);
1356 :
1357 : // add shape fct at ip * value at shape
1358 0 : approxSolIP += valSH * rTrialSpace.shape(sh, vLocIP[ip]);
1359 :
1360 : // add gradient at ip
1361 : VecScaleAppend(locTmp, valSH, vLocGradient[sh]);
1362 : }
1363 :
1364 : // compute global gradient
1365 : MathVector<worldDim> approxGradIP;
1366 : MathMatrix<worldDim, elemDim> JTInv;
1367 0 : map.jacobian_transposed_inverse(JTInv, vLocIP[ip]); //Inverse(JTInv, vJT[ip]);
1368 :
1369 : #ifdef UG_DEBUG
1370 : CheckGeneralizedInverse<worldDim,elemDim>(vJT[ip], JTInv);
1371 : #endif //UG_DEBUG
1372 :
1373 : MatVecMult(approxGradIP, JTInv, locTmp);
1374 :
1375 : // get squared of difference
1376 0 : vValue[ip] = (exactSolIP - approxSolIP) * (exactSolIP - approxSolIP);
1377 0 : vValue[ip] += VecDistanceSq(approxGradIP, exactGradIP);
1378 : }
1379 :
1380 0 : }
1381 0 : UG_CATCH_THROW("H1ErrorIntegrand::evaluate: trial space missing.");
1382 0 : }
1383 : };
1384 :
1385 : /// compute H1 error of a function on the whole domain or on some subsets
1386 : /**
1387 : * This function computes the H1-difference between a given analytic function
1388 : * and a grid function.
1389 : *
1390 : * \param[in] spExactSol analytic function
1391 : * \param[in] spExactGrad analytic gradient
1392 : * \param[in] spGridFct grid function
1393 : * \param[in] cmp symbolic name of component function
1394 : * \param[in] time time point
1395 : * \param[in] quadOrder order of quadrature rule
1396 : * \param[in] subsets subsets, where to compute
1397 : * \returns number l2-norm of difference
1398 : */
1399 : template <typename TGridFunction>
1400 0 : number H1Error(SmartPtr<UserData<number, TGridFunction::dim> > spExactSol,
1401 : SmartPtr<UserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > spExactGrad,
1402 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
1403 : number time, int quadOrder, const char* subsets)
1404 : {
1405 : // get function id of name
1406 : const size_t fct = spGridFct->fct_id_by_name(cmp);
1407 :
1408 : // check that function exists
1409 0 : if(fct >= spGridFct->num_fct())
1410 0 : UG_THROW("H1Error: Function space does not contain"
1411 : " a function with name " << cmp << ".");
1412 :
1413 0 : H1ErrorIntegrand<TGridFunction> spIntegrand(spExactSol, spExactGrad, *spGridFct, fct, time);
1414 0 : return sqrt(IntegrateSubsets(spIntegrand, *spGridFct, subsets, quadOrder));
1415 : }
1416 :
1417 : template <typename TGridFunction>
1418 0 : number H1Error(SmartPtr<UserData<number, TGridFunction::dim> > spExactSol,
1419 : SmartPtr<UserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > spExactGrad,
1420 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
1421 : number time, int quadOrder)
1422 : {
1423 0 : return H1Error(spExactSol, spExactGrad, spGridFct, cmp, time, quadOrder, NULL);
1424 : }
1425 :
1426 : #ifdef UG_FOR_LUA
1427 : template <typename TGridFunction>
1428 0 : number H1Error(const char* ExactSol, const char* ExactGrad,
1429 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
1430 : number time, int quadOrder, const char* subsets)
1431 : {
1432 : static const int dim = TGridFunction::domain_type::dim;
1433 0 : SmartPtr<UserData<number, dim> > spExactSol
1434 0 : = make_sp(new LuaUserData<number, dim>(ExactSol));
1435 0 : SmartPtr<UserData<MathVector<dim>, dim> > spExactGrad
1436 0 : = make_sp(new LuaUserData<MathVector<dim>, dim>(ExactGrad));
1437 0 : return H1Error(spExactSol, spExactGrad, spGridFct, cmp, time, quadOrder, subsets);
1438 : }
1439 :
1440 : template <typename TGridFunction>
1441 0 : number H1Error(const char* ExactSol, const char* ExactGrad,
1442 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
1443 : number time, int quadOrder)
1444 : {
1445 0 : return H1Error(ExactSol, ExactGrad, spGridFct, cmp, time, quadOrder, NULL);
1446 : }
1447 : #endif
1448 :
1449 :
1450 :
1451 : /// computes an (abstract) distance between two functions
1452 : /**
1453 : * This function computes the (abstract) TDiffError-difference between two grid functions that
1454 : * may be defined on different grids. The element loop is performed over the
1455 : * finer level.
1456 : *
1457 : * \param[in] spGridFct1 grid function 1
1458 : * \param[in] cmp1 symbolic name of component function
1459 : * \param[in] spGridFct2 grid function 2
1460 : * \param[in] cmp2 symbolic name of component function
1461 : * \param[in] quadOrder order of quadrature rule
1462 : * \param[in] subsets subsets, where to compute
1463 : * \returns number H1-norm of difference
1464 : */
1465 : template <typename TDistIntegrand, typename TGridFunction>
1466 0 : number GridFunctionDistance2(TGridFunction& spGridFct1, const char* cmp1,
1467 : TGridFunction& spGridFct2, const char* cmp2,
1468 : int quadOrder, const char* subsets)
1469 : {
1470 : // get function id of name
1471 : const size_t fct1 = spGridFct1.fct_id_by_name(cmp1);
1472 : const size_t fct2 = spGridFct2.fct_id_by_name(cmp2);
1473 :
1474 : // check that function exists
1475 0 : if(fct1 >= spGridFct1.num_fct())
1476 0 : UG_THROW("GridFunctionDistance: Function space does not contain"
1477 : " a function with name " << cmp1 << ".");
1478 0 : if(fct2 >= spGridFct2.num_fct())
1479 0 : UG_THROW("GridFunctionDistance: Function space does not contain"
1480 : " a function with name " << cmp2 << ".");
1481 :
1482 : // get top level of grid functions
1483 0 : const int level1 = spGridFct1.dof_distribution()->grid_level().level();
1484 0 : const int level2 = spGridFct2.dof_distribution()->grid_level().level();
1485 :
1486 0 : if(level1 > level2){
1487 : // level check
1488 0 : TDistIntegrand spIntegrand(spGridFct1, fct1, spGridFct2, fct2);
1489 0 : return IntegrateSubsets(spIntegrand, spGridFct1, subsets, quadOrder);
1490 : }else{
1491 0 : TDistIntegrand spIntegrand(spGridFct2, fct2, spGridFct1, fct1);
1492 0 : return IntegrateSubsets(spIntegrand, spGridFct2, subsets, quadOrder);
1493 : }
1494 :
1495 : }
1496 :
1497 : //! Computes (weighted) distance
1498 : template <typename TDistIntegrand, typename TGridFunction>
1499 0 : number GridFunctionDistance2(TGridFunction& spGridFct1, const char* cmp1,
1500 : TGridFunction& spGridFct2, const char* cmp2,
1501 : int quadOrder, const char* subsets, ConstSmartPtr<typename TDistIntegrand::weight_type> spWeights)
1502 : {
1503 : // get function id of name
1504 : const size_t fct1 = spGridFct1.fct_id_by_name(cmp1);
1505 : const size_t fct2 = spGridFct2.fct_id_by_name(cmp2);
1506 :
1507 : // check that function exists
1508 0 : if(fct1 >= spGridFct1.num_fct())
1509 0 : UG_THROW("GridFunctionDistance: Function space does not contain"
1510 : " a function with name " << cmp1 << ".");
1511 0 : if(fct2 >= spGridFct2.num_fct())
1512 0 : UG_THROW("GridFunctionDistance: Function space does not contain"
1513 : " a function with name " << cmp2 << ".");
1514 :
1515 : // get top level of gridfunctions
1516 0 : const int level1 = spGridFct1.dof_distribution()->grid_level().level();
1517 0 : const int level2 = spGridFct2.dof_distribution()->grid_level().level();
1518 :
1519 : // w/ weights
1520 0 : if(level1 > level2){
1521 0 : TDistIntegrand spIntegrand(spGridFct1, fct1, spGridFct2, fct2, spWeights);
1522 0 : return IntegrateSubsets(spIntegrand, spGridFct1, subsets, quadOrder);
1523 : }else{
1524 0 : TDistIntegrand spIntegrand(spGridFct2, fct2, spGridFct1, fct1, spWeights);
1525 0 : return IntegrateSubsets(spIntegrand, spGridFct2, subsets, quadOrder);
1526 : }
1527 : }
1528 :
1529 :
1530 : //! Computes (weighted) distance with shift for averages
1531 : template <typename TDistIntegrand, typename TGridFunction>
1532 0 : number GridFunctionDistance2(TGridFunction& spGridFct1, const char* cmp1,
1533 : TGridFunction& spGridFct2, const char* cmp2,
1534 : int quadOrder, const char* subsets,
1535 : ConstSmartPtr<typename TDistIntegrand::weight_type> spWeights,
1536 : number distAvg12)
1537 : {
1538 : // get function id of name
1539 : const size_t fct1 = spGridFct1.fct_id_by_name(cmp1);
1540 : const size_t fct2 = spGridFct2.fct_id_by_name(cmp2);
1541 :
1542 : // check that function exists
1543 0 : if(fct1 >= spGridFct1.num_fct())
1544 0 : UG_THROW("GridFunctionDistance: Function space does not contain"
1545 : " a function with name " << cmp1 << ".");
1546 0 : if(fct2 >= spGridFct2.num_fct())
1547 0 : UG_THROW("GridFunctionDistance: Function space does not contain"
1548 : " a function with name " << cmp2 << ".");
1549 :
1550 : // get top level of gridfunctions
1551 0 : const int level1 = spGridFct1.dof_distribution()->grid_level().level();
1552 0 : const int level2 = spGridFct2.dof_distribution()->grid_level().level();
1553 :
1554 : // w/ weights
1555 0 : if(level1 > level2){
1556 0 : TDistIntegrand spIntegrand(spGridFct1, fct1, spGridFct2, fct2, spWeights, distAvg12);
1557 0 : return IntegrateSubsets(spIntegrand, spGridFct1, subsets, quadOrder);
1558 : }else{
1559 0 : TDistIntegrand spIntegrand(spGridFct2, fct2, spGridFct1, fct1, spWeights, -distAvg12);
1560 0 : return IntegrateSubsets(spIntegrand, spGridFct2, subsets, quadOrder);
1561 : }
1562 : }
1563 :
1564 :
1565 :
1566 : ////////////////////////////////////////////////////////////////////////////////
1567 : // L2 Integrand
1568 : ////////////////////////////////////////////////////////////////////////////////
1569 :
1570 : /// Grid function as L2 integrand
1571 : template <typename TGridFunction>
1572 : class L2Integrand
1573 : : public StdIntegrand<number, TGridFunction::dim, L2Integrand<TGridFunction> >
1574 : {
1575 : public:
1576 : // world dimension of grid function
1577 : static const int worldDim = TGridFunction::dim;
1578 : typedef UserData<number, worldDim> weight_type;
1579 :
1580 : protected:
1581 : // grid function data
1582 : ScalarGridFunctionData<TGridFunction> m_scalarData;
1583 :
1584 : /// scalar weight (optional, default is 1.0)
1585 : ConstSmartPtr<weight_type> m_spWeight;
1586 :
1587 : public:
1588 : /// CTOR
1589 0 : L2Integrand(TGridFunction& spGridFct, size_t cmp)
1590 0 : : m_scalarData(spGridFct, cmp), m_spWeight(make_sp(new ConstUserNumber<worldDim>(1.0)))
1591 0 : {};
1592 :
1593 0 : L2Integrand(TGridFunction& spGridFct, size_t cmp, ConstSmartPtr<weight_type> spWeight)
1594 0 : : m_scalarData(spGridFct, cmp), m_spWeight(spWeight)
1595 : {};
1596 :
1597 : /// DTOR
1598 0 : virtual ~L2Integrand() {};
1599 :
1600 : /// sets subset
1601 0 : virtual void set_subset(int si)
1602 : {
1603 0 : if(!m_scalarData.is_def_in_subset(si))
1604 0 : UG_THROW("L2ErrorIntegrand: Grid function component" <<m_scalarData.fct() <<" not defined on subset "<<si);
1605 : IIntegrand<number, worldDim>::set_subset(si);
1606 0 : }
1607 :
1608 : /// \copydoc IIntegrand::values
1609 : template <int elemDim>
1610 0 : void evaluate(number vValue[],
1611 : const MathVector<worldDim> vGlobIP[],
1612 : GridObject* pElem,
1613 : const MathVector<worldDim> vCornerCoords[],
1614 : const MathVector<elemDim> vLocIP[],
1615 : const MathMatrix<elemDim, worldDim> vJT[],
1616 : const size_t numIP)
1617 : {
1618 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
1619 0 : ReferenceObjectID roid = (ReferenceObjectID) pElem->reference_object_id();
1620 :
1621 : // element weights
1622 : typedef typename weight_type::data_type ipdata_type;
1623 0 : std::vector<ipdata_type> locElemWeights(numIP, 1.0);
1624 : UG_ASSERT(m_spWeight.valid(), "L2Integrand::evaluate requires valid weights!");
1625 0 : (*m_spWeight)(&locElemWeights[0], vGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), numIP);
1626 :
1627 : try{
1628 : // get trial space
1629 : const LocalShapeFunctionSet<elemDim>& rTrialSpace =
1630 : LocalFiniteElementProvider::get<elemDim>(roid, m_scalarData.id());
1631 :
1632 : // number of dofs on element
1633 0 : const size_t num_sh = rTrialSpace.num_sh();
1634 :
1635 : // get multiindices of element
1636 : std::vector<DoFIndex> ind; // aux. index array
1637 : m_scalarData.dof_indices(pElem, ind);
1638 :
1639 : // check multi indices
1640 0 : if(ind.size() != num_sh)
1641 0 : UG_THROW("L2Integrand::evaluate: Wrong number of multi indices.");
1642 :
1643 : // loop all integration points
1644 0 : for(size_t ip = 0; ip < numIP; ++ip)
1645 : {
1646 :
1647 : // compute approximated solution at integration point
1648 : number approxSolIP = 0.0;
1649 0 : for(size_t sh = 0; sh < num_sh; ++sh)
1650 : {
1651 : // get value at shape point (e.g. corner for P1 fct)
1652 : // and add shape fct at ip * value at shape
1653 0 : const number valSH = DoFRef(m_scalarData.grid_function(), ind[sh]);
1654 0 : approxSolIP += valSH * rTrialSpace.shape(sh, vLocIP[ip]);
1655 : }
1656 :
1657 : // get square
1658 0 : vValue[ip] = locElemWeights[ip]*approxSolIP*approxSolIP;
1659 :
1660 : }
1661 :
1662 : }
1663 0 : UG_CATCH_THROW("L2FuncIntegrand::values: trial space missing.");
1664 0 : }
1665 : };
1666 :
1667 : /**
1668 : * This function computes the square of the L2-norm of a grid function.
1669 : *
1670 : * \param[in] spGridFct grid function
1671 : * \param[in] cmp symbolic name of function
1672 : * \param[in] quadOrder order of quadrature rule
1673 : * \param[in] subsets subsets, where to interpolate
1674 : * (NULL indicates that all full-dimensional subsets
1675 : * shall be considered)
1676 : * \returns number l2-norm
1677 : */
1678 :
1679 : template <typename TGridFunction>
1680 0 : number L2Norm2(TGridFunction& u, const char* cmp,
1681 : int quadOrder, const char* subsets,
1682 : ConstSmartPtr<typename L2Integrand<TGridFunction>::weight_type> spWeight)
1683 : {
1684 : // get function id of name
1685 : const size_t fct = u.fct_id_by_name(cmp);
1686 :
1687 : // check that function exists
1688 0 : UG_COND_THROW(fct >= u.num_fct(), "L2Norm: Function space does not contain"
1689 : " a function with name " << cmp << ".");
1690 :
1691 0 : L2Integrand<TGridFunction> integrandL2(u, fct, spWeight);
1692 0 : return IntegrateSubsets(integrandL2, u, subsets, quadOrder);
1693 : }
1694 :
1695 : template <typename TGridFunction>
1696 0 : number L2Norm2(TGridFunction& u, const char* cmp,
1697 : int quadOrder, const char* subsets)
1698 : {
1699 : // get function id of name
1700 : const size_t fct = u.fct_id_by_name(cmp);
1701 :
1702 : // check that function exists
1703 0 : UG_COND_THROW(fct >= u.num_fct(), "L2Norm: Function space does not contain"
1704 : " a function with name " << cmp << ".");
1705 :
1706 0 : L2Integrand<TGridFunction> integrandL2(u, fct);
1707 0 : return IntegrateSubsets(integrandL2, u, subsets, quadOrder);
1708 : }
1709 :
1710 : template <typename TGridFunction>
1711 : number L2Norm(TGridFunction& u, const char* cmp,
1712 : int quadOrder, const char* subsets)
1713 : {
1714 0 : return sqrt(L2Norm2(u, cmp, quadOrder, subsets));
1715 : }
1716 : /**
1717 : * This function computes the L2-norm of a grid function on all full-dim subsets.
1718 : *
1719 : * \param[in] spGridFct grid function
1720 : * \param[in] cmp symbolic name of function
1721 : * \param[in] quadOrder order of quadrature rule
1722 : * \returns number l2-norm
1723 : */
1724 : template <typename TGridFunction>
1725 : number L2Norm(TGridFunction& gridFct, const char* cmp, int quadOrder)
1726 : { return L2Norm(gridFct, cmp, quadOrder, NULL); }
1727 :
1728 : template <typename TGridFunction>
1729 0 : number L2Norm(SmartPtr<TGridFunction> spGridFct, const char* cmp, int quadOrder, const char* subsets)
1730 0 : { return L2Norm(*spGridFct, cmp, quadOrder, subsets); }
1731 :
1732 : template <typename TGridFunction>
1733 0 : number L2Norm(SmartPtr<TGridFunction> spGridFct, const char* cmp, int quadOrder)
1734 0 : { return L2Norm(spGridFct, cmp, quadOrder, NULL); }
1735 :
1736 :
1737 : /// Integrand for the distance of two grid functions - evaluated in the (weighted) H1-semi norm
1738 : template <typename TGridFunction>
1739 : class L2DistIntegrand
1740 : : public StdIntegrand<number, TGridFunction::dim, L2DistIntegrand<TGridFunction> >
1741 : {
1742 : public:
1743 : /// world dimension of grid function
1744 : static const int worldDim = TGridFunction::dim;
1745 : typedef typename L2Integrand<TGridFunction>::weight_type weight_type;
1746 :
1747 : protected:
1748 : ScalarGridFunctionData<TGridFunction> m_fineData;
1749 : const int m_fineTopLevel;
1750 :
1751 : ScalarGridFunctionData<TGridFunction> m_coarseData;
1752 : const int m_coarseTopLevel;
1753 :
1754 : /// multigrid
1755 : SmartPtr<MultiGrid> m_spMG;
1756 :
1757 : ConstSmartPtr<weight_type> m_spWeight;
1758 :
1759 : /// shift
1760 : double m_deltaFineCoarse;
1761 :
1762 : public:
1763 :
1764 : /// constructor (1st is fine grid function)
1765 0 : L2DistIntegrand(TGridFunction& fineGridFct, size_t fineCmp,
1766 : TGridFunction& coarseGridFct, size_t coarseCmp)
1767 0 : : m_fineData(fineGridFct, fineCmp), m_fineTopLevel(fineGridFct.dof_distribution()->grid_level().level()),
1768 0 : m_coarseData(coarseGridFct, coarseCmp), m_coarseTopLevel(coarseGridFct.dof_distribution()->grid_level().level()),
1769 0 : m_spMG(m_fineData.domain()->grid()), m_spWeight(make_sp(new ConstUserNumber<TGridFunction::dim>(1.0))),
1770 0 : m_deltaFineCoarse(0.0)
1771 : {
1772 0 : UG_COND_THROW(m_fineTopLevel < m_coarseTopLevel,
1773 : "L2DiffIntegrand: fine and top level inverted.");
1774 0 : UG_COND_THROW(m_fineData.domain().get() != m_coarseData.domain().get(),
1775 : "L2DiffIntegrand: grid functions defined on different domains.");
1776 0 : };
1777 :
1778 : /// constructor (1st is fine grid function)
1779 : L2DistIntegrand(TGridFunction& fineGridFct, size_t fineCmp,
1780 : TGridFunction& coarseGridFct, size_t coarseCmp, ConstSmartPtr<weight_type> spWeight)
1781 : : m_fineData(fineGridFct, fineCmp), m_fineTopLevel(fineGridFct.dof_distribution()->grid_level().level()),
1782 : m_coarseData(coarseGridFct, coarseCmp), m_coarseTopLevel(coarseGridFct.dof_distribution()->grid_level().level()),
1783 : m_spMG(m_fineData.domain()->grid()), m_spWeight(spWeight),
1784 : m_deltaFineCoarse(0.0)
1785 : {
1786 : UG_COND_THROW(m_fineTopLevel < m_coarseTopLevel,
1787 : "L2DiffIntegrand: fine and top level inverted.");
1788 : UG_COND_THROW(m_fineData.domain().get() != m_coarseData.domain().get(),
1789 : "L2DiffIntegrand: grid functions defined on different domains.");
1790 : };
1791 :
1792 : /// constructor (1st is fine grid function)
1793 0 : L2DistIntegrand(TGridFunction& fineGridFct, size_t fineCmp,
1794 : TGridFunction& coarseGridFct, size_t coarseCmp, ConstSmartPtr<weight_type> spWeight, number dist12)
1795 0 : : m_fineData(fineGridFct, fineCmp), m_fineTopLevel(fineGridFct.dof_distribution()->grid_level().level()),
1796 0 : m_coarseData(coarseGridFct, coarseCmp), m_coarseTopLevel(coarseGridFct.dof_distribution()->grid_level().level()),
1797 0 : m_spMG(m_fineData.domain()->grid()), m_spWeight(spWeight),
1798 0 : m_deltaFineCoarse(dist12)
1799 : {
1800 0 : UG_COND_THROW(m_fineTopLevel < m_coarseTopLevel,
1801 : "L2DiffIntegrand: fine and top level inverted.");
1802 0 : UG_COND_THROW(m_fineData.domain().get() != m_coarseData.domain().get(),
1803 : "L2DiffIntegrand: grid functions defined on different domains.");
1804 0 : };
1805 :
1806 :
1807 0 : virtual ~L2DistIntegrand() {}
1808 :
1809 : /// sets subset
1810 0 : virtual void set_subset(int si)
1811 : {
1812 0 : UG_COND_THROW(!m_fineData.is_def_in_subset(si),
1813 : "L2DiffIntegrand: Grid function component" <<m_fineData.fct()<<" not defined on subset "<<si);
1814 0 : UG_COND_THROW(!m_coarseData.is_def_in_subset(si),
1815 : "L2DiffIntegrand: Grid function component" <<m_coarseData.fct()<<" not defined on subset "<<si);
1816 : IIntegrand<number, worldDim>::set_subset(si);
1817 0 : }
1818 :
1819 : /// \copydoc IIntegrand::values
1820 : template <int elemDim>
1821 0 : void evaluate(number vValue[],
1822 : const MathVector<worldDim> vFineGlobIP[],
1823 : GridObject* pFineElem,
1824 : const MathVector<worldDim> vCornerCoords[],
1825 : const MathVector<elemDim> vFineLocIP[],
1826 : const MathMatrix<elemDim, worldDim> vJT[],
1827 : const size_t numIP)
1828 : {
1829 : typedef typename TGridFunction::template dim_traits<elemDim>::grid_base_object Element;
1830 :
1831 : // get coarse element
1832 : GridObject* pCoarseElem = pFineElem;
1833 0 : if(m_coarseTopLevel < m_fineTopLevel){
1834 : int parentLevel = m_spMG->get_level(pCoarseElem);
1835 0 : while(parentLevel > m_coarseTopLevel){
1836 0 : pCoarseElem = m_spMG->get_parent(pCoarseElem);
1837 : parentLevel = m_spMG->get_level(pCoarseElem);
1838 : }
1839 : }
1840 :
1841 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
1842 0 : const ReferenceObjectID fineROID = pFineElem->reference_object_id();
1843 0 : const ReferenceObjectID coarseROID = pCoarseElem->reference_object_id();
1844 :
1845 : // get corner coordinates
1846 : std::vector<MathVector<worldDim> > vCornerCoarse;
1847 0 : CollectCornerCoordinates(vCornerCoarse, *static_cast<Element*>(pCoarseElem), *m_coarseData.domain());
1848 :
1849 : // get Reference Mapping
1850 : DimReferenceMapping<elemDim, worldDim>& map
1851 : = ReferenceMappingProvider::get<elemDim, worldDim>(coarseROID, vCornerCoarse);
1852 :
1853 : std::vector<MathVector<elemDim> > vCoarseLocIP;
1854 0 : vCoarseLocIP.resize(numIP);
1855 0 : for(size_t ip = 0; ip < vCoarseLocIP.size(); ++ip) VecSet(vCoarseLocIP[ip], 0.0);
1856 0 : map.global_to_local(&vCoarseLocIP[0], vFineGlobIP, numIP);
1857 :
1858 : // element weights
1859 : typedef typename weight_type::data_type ipdata_type;
1860 0 : std::vector<ipdata_type> fineElemWeights(numIP, 1.0);
1861 : UG_ASSERT(m_spWeight.valid(), "L2DistIntegrand::evaluate requires valid weights! ");
1862 0 : (*m_spWeight)(&fineElemWeights[0], vFineGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), numIP);
1863 :
1864 : try{
1865 : // get trial space
1866 : const LocalShapeFunctionSet<elemDim>& rFineLSFS =
1867 : LocalFiniteElementProvider::get<elemDim>(fineROID, m_fineData.id());
1868 : const LocalShapeFunctionSet<elemDim>& rCoarseLSFS =
1869 : LocalFiniteElementProvider::get<elemDim>(coarseROID, m_coarseData.id());
1870 :
1871 : // get multiindices of element
1872 : std::vector<DoFIndex> vFineMI, vCoarseMI;
1873 : m_fineData.dof_indices(pFineElem, vFineMI);
1874 : m_coarseData.dof_indices(pCoarseElem, vCoarseMI);
1875 :
1876 : // loop all integration points
1877 0 : for(size_t ip = 0; ip < numIP; ++ip)
1878 : {
1879 : // compute approximated solution at integration point
1880 : number fineSolIP = 0.0;
1881 0 : for(size_t sh = 0; sh < vFineMI.size(); ++sh)
1882 : {
1883 : // get value at shape point (e.g. corner for P1 fct)
1884 0 : const number val = DoFRef(m_fineData.grid_function(), vFineMI[sh]);
1885 :
1886 : // add shape fct at ip * value at shape
1887 0 : fineSolIP += val * rFineLSFS.shape(sh, vFineLocIP[ip]);
1888 : }
1889 : number coarseSolIP = 0.0;
1890 0 : for(size_t sh = 0; sh < vCoarseMI.size(); ++sh)
1891 : {
1892 : // get value at shape point (e.g. corner for P1 fct)
1893 0 : const number val = DoFRef(m_coarseData.grid_function(), vCoarseMI[sh]);
1894 :
1895 : // add shape fct at ip * value at shape
1896 0 : coarseSolIP += val * rCoarseLSFS.shape(sh, vCoarseLocIP[ip]);
1897 : }
1898 :
1899 : // get squared of difference
1900 0 : vValue[ip] = fineElemWeights[ip]*(fineSolIP - coarseSolIP -m_deltaFineCoarse)*(fineSolIP-coarseSolIP-m_deltaFineCoarse);
1901 : }
1902 :
1903 : }
1904 0 : UG_CATCH_THROW("L2DistIntegrand::evaluate: trial space missing.");
1905 0 : }
1906 : };
1907 :
1908 :
1909 : /// computes the squared l2 distance between two functions
1910 : template <typename TGridFunction>
1911 0 : number L2Distance2(TGridFunction& spGridFct1, const char* cmp1,
1912 : TGridFunction& spGridFct2, const char* cmp2,
1913 : int quadOrder, const char* subsets,
1914 : ConstSmartPtr<typename L2Integrand<TGridFunction>::weight_type> spWeight, number avgDist12=0.0)
1915 : {
1916 : return GridFunctionDistance2<L2DistIntegrand<TGridFunction>, TGridFunction>
1917 0 : (spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, subsets, spWeight, avgDist12);
1918 : }
1919 :
1920 :
1921 : /// computes the squared l2 distance between two functions
1922 : template <typename TGridFunction>
1923 : number L2Distance2(TGridFunction& spGridFct1, const char* cmp1,
1924 : TGridFunction& spGridFct2, const char* cmp2,
1925 : int quadOrder, const char* subsets)
1926 : {
1927 : return GridFunctionDistance2<L2DistIntegrand<TGridFunction>, TGridFunction>
1928 0 : (spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, subsets);
1929 : }
1930 :
1931 : /// computes the l2 distance between two functions
1932 : template <typename TGridFunction>
1933 : number L2Distance(TGridFunction& spGridFct1, const char* cmp1,
1934 : TGridFunction& spGridFct2, const char* cmp2,
1935 : int quadOrder, const char* subsets)
1936 : {
1937 0 : return sqrt(L2Distance2(spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, subsets));
1938 : }
1939 :
1940 :
1941 : template <typename TGridFunction>
1942 0 : number L2Error(SmartPtr<TGridFunction> spGridFct1, const char* cmp1,
1943 : SmartPtr<TGridFunction> spGridFct2, const char* cmp2,
1944 : int quadOrder)
1945 : {
1946 0 : return L2Distance(*spGridFct1, cmp1, *spGridFct2, cmp2, quadOrder, NULL);
1947 : }
1948 :
1949 : template <typename TGridFunction>
1950 0 : number L2Error(SmartPtr<TGridFunction> spGridFct1, const char* cmp1,
1951 : SmartPtr<TGridFunction> spGridFct2, const char* cmp2,
1952 : int quadOrder, const char* subsets)
1953 : {
1954 0 : return L2Distance(*spGridFct1, cmp1, *spGridFct2, cmp2, quadOrder, subsets);
1955 : }
1956 :
1957 : ////////////////////////////////////////////////////////////////////////////////
1958 : // H1 semi-norm Integrand
1959 : ////////////////////////////////////////////////////////////////////////////////
1960 :
1961 :
1962 :
1963 : //! Norm of a grid function, evaluated in (weighted) H1-semi norm
1964 : template <typename TGridFunction>
1965 : class H1SemiIntegrand
1966 : : public StdIntegrand<number, TGridFunction::dim, H1SemiIntegrand<TGridFunction> >
1967 : {
1968 : public:
1969 : /// world dimension of grid function
1970 : static const int worldDim = TGridFunction::dim;
1971 : typedef UserData<MathMatrix<worldDim, worldDim>, worldDim> weight_type;
1972 :
1973 : protected:
1974 : /// grid function data
1975 : ScalarGridFunctionData<TGridFunction> m_scalarData;
1976 :
1977 : /// scalar weight (optional)
1978 : ConstSmartPtr<weight_type> m_spWeight;
1979 :
1980 : public:
1981 : /// constructor
1982 0 : H1SemiIntegrand(TGridFunction& gridFct, size_t cmp)
1983 0 : : m_scalarData(gridFct, cmp), m_spWeight(make_sp(new ConstUserMatrix<worldDim>(1.0))) {}
1984 :
1985 : /// constructor
1986 0 : H1SemiIntegrand(TGridFunction& gridFct, size_t cmp, ConstSmartPtr<weight_type> spWeight)
1987 0 : : m_scalarData(gridFct, cmp), m_spWeight(spWeight) {}
1988 :
1989 : /// DTOR
1990 0 : virtual ~H1SemiIntegrand() {};
1991 :
1992 : /// sets subset
1993 0 : virtual void set_subset(int si)
1994 : {
1995 :
1996 0 : UG_COND_THROW(!m_scalarData.is_def_in_subset(si), "H1Error: Grid function component"
1997 : <<m_scalarData.fct()<<" not defined on subset "<<si);
1998 : IIntegrand<number, worldDim>::set_subset(si);
1999 0 : }
2000 :
2001 : /// \copydoc IIntegrand::values
2002 : template <int elemDim>
2003 0 : void evaluate(number vValue[],
2004 : const MathVector<worldDim> vGlobIP[],
2005 : GridObject* pElem,
2006 : const MathVector<worldDim> vCornerCoords[],
2007 : const MathVector<elemDim> vLocIP[],
2008 : const MathMatrix<elemDim, worldDim> vJT[],
2009 : const size_t numIP)
2010 : {
2011 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
2012 0 : const ReferenceObjectID roid = pElem->reference_object_id();
2013 : const TGridFunction &gridFct= m_scalarData.grid_function();
2014 :
2015 : DimReferenceMapping<elemDim, worldDim>& map
2016 : = ReferenceMappingProvider::get<elemDim, worldDim>(roid, vCornerCoords);
2017 :
2018 : typedef typename weight_type::data_type ipdata_type;
2019 :
2020 0 : std::vector<ipdata_type> elemWeights(numIP, MathMatrix<worldDim, worldDim>());
2021 : UG_ASSERT(m_spWeight.valid(), "H1SemiIntegrand::evaluate requires valid weights!");
2022 :
2023 :
2024 0 : if(m_spWeight->requires_grid_fct())
2025 : {
2026 : // get local solution if needed
2027 : LocalIndices ind;
2028 : LocalVector uloc;
2029 : gridFct.indices(pElem, ind); // get global indices
2030 0 : uloc.resize(ind); // adapt local algebra
2031 0 : GetLocalVector(uloc, gridFct); // read local values of u
2032 :
2033 : // compute data
2034 : try{
2035 0 : (*m_spWeight)(&elemWeights[0], vGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), pElem,
2036 : vCornerCoords, vLocIP, numIP, &uloc, NULL);
2037 0 : } UG_CATCH_THROW("H1SemiIntegrand: Cannot evaluate weight data.");
2038 : }
2039 : else
2040 : {
2041 : // compute data
2042 : try{
2043 0 : (*m_spWeight)(&elemWeights[0], vGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), numIP);
2044 0 : } UG_CATCH_THROW("H1SemiIntegrand: Cannot evaluate weight data.");
2045 : }
2046 :
2047 :
2048 : try{
2049 :
2050 :
2051 :
2052 :
2053 : // get trial space
2054 : const LocalShapeFunctionSet<elemDim>& rTrialSpace =
2055 : LocalFiniteElementProvider::get<elemDim>(roid, m_scalarData.id());
2056 :
2057 : // number of dofs on element
2058 0 : const size_t num_sh = rTrialSpace.num_sh();
2059 :
2060 : // get multiindices of element
2061 : std::vector<DoFIndex> ind; // aux. index array
2062 : gridFct.dof_indices(pElem, m_scalarData.fct(), ind);
2063 :
2064 : // check multi indices
2065 0 : UG_COND_THROW(ind.size() != num_sh, "H1SemiNormFuncIntegrand::evaluate: Wrong number of multi-)indices.");
2066 :
2067 : // loop all integration points
2068 0 : std::vector<MathVector<elemDim> > vLocGradient(num_sh);
2069 0 : for(size_t ip = 0; ip < numIP; ++ip)
2070 : {
2071 : // compute shape gradients at ip
2072 0 : rTrialSpace.grads(&vLocGradient[0], vLocIP[ip]);
2073 :
2074 : // compute approximated solution at integration point
2075 : number approxSolIP = 0.0;
2076 : MathVector<elemDim> tmpVec(0.0);
2077 0 : for(size_t sh = 0; sh < num_sh; ++sh)
2078 : {
2079 : // get value at shape point (e.g. corner for P1 fct)
2080 0 : const number valSH = DoFRef(gridFct, ind[sh]);
2081 :
2082 : // add shape fct at ip * value at shape
2083 0 : approxSolIP += valSH * rTrialSpace.shape(sh, vLocIP[ip]);
2084 :
2085 : // add gradient at ip
2086 : VecScaleAppend(tmpVec, valSH, vLocGradient[sh]);
2087 : }
2088 :
2089 : // compute gradient
2090 : MathVector<worldDim> approxGradIP;
2091 : MathMatrix<worldDim, elemDim> JTInv;
2092 0 : map.jacobian_transposed_inverse(JTInv, vLocIP[ip]);//Inverse(JTInv, vJT[ip]);
2093 :
2094 : #ifdef UG_DEBUG
2095 : CheckGeneralizedInverse<worldDim,elemDim>(vJT[ip], JTInv);
2096 : #endif //UG_DEBUG
2097 :
2098 : MatVecMult(approxGradIP, JTInv, tmpVec);
2099 :
2100 : // get norm squared
2101 : MathVector<worldDim> approxDGradIP;
2102 : MatVecMult(approxDGradIP, elemWeights[ip], approxGradIP);
2103 0 : vValue[ip] = VecDot(approxDGradIP, approxGradIP);
2104 : }
2105 :
2106 0 : }
2107 0 : UG_CATCH_THROW("H1SemiIntegrand::evaluate: trial space missing.");
2108 0 : }
2109 : };
2110 :
2111 :
2112 : /// compute H1 semi-norm of a function on the whole domain (or on selected subsets)
2113 : /**
2114 : * This function computes the H1 semi-norm of a grid function
2115 : *
2116 : * \param[in] spGridFct grid function
2117 : * \param[in] cmp symbolic name of component function
2118 : * \param[in] time time point
2119 : * \param[in] quadOrder order of quadrature rule
2120 : * \param[in] subsets subsets, where to compute (OPTIONAL)
2121 : * \param[in] weights element-wise weights (OPTIONAL)
2122 : * \returns number l2-norm
2123 : */
2124 : template <typename TGridFunction>
2125 0 : number H1SemiNorm2(TGridFunction& gridFct, const char* cmp, int quadOrder, const char* subsets=NULL,
2126 : ConstSmartPtr<typename H1SemiIntegrand<TGridFunction>::weight_type> weights = SPNULL)
2127 : {
2128 : // get function id of name
2129 : const size_t fct = gridFct.fct_id_by_name(cmp);
2130 :
2131 : // check that function exists
2132 0 : UG_COND_THROW(fct >= gridFct.num_fct(), "H1SemiNorm: Function space does not contain"
2133 : " a function with name " << cmp << ".");
2134 0 : if (weights.invalid()) {
2135 0 : H1SemiIntegrand<TGridFunction> integrand(gridFct, fct);
2136 0 : return IntegrateSubsets(integrand, gridFct, subsets, quadOrder);
2137 : } else {
2138 0 : H1SemiIntegrand<TGridFunction> integrand(gridFct, fct, weights);
2139 0 : return IntegrateSubsets(integrand, gridFct, subsets, quadOrder);
2140 : }
2141 : }
2142 :
2143 : /// Computes the H1SemiNorm
2144 : template <typename TGridFunction>
2145 0 : number H1SemiNorm(TGridFunction& gridFct, const char* cmp, int quadOrder, const char* subsets=NULL,
2146 : ConstSmartPtr<typename H1SemiIntegrand<TGridFunction>::weight_type> weights = SPNULL)
2147 : {
2148 0 : return (sqrt(H1SemiNorm2(gridFct, cmp, quadOrder, subsets, weights)));
2149 : }
2150 :
2151 : // Delegating to H1SemiNorm
2152 : template <typename TGridFunction>
2153 0 : number H1SemiNorm(SmartPtr<TGridFunction> spGridFct, const char* cmp, int quadOrder, const char* subsets)
2154 0 : { return H1SemiNorm(*spGridFct, cmp, quadOrder, subsets); }
2155 :
2156 : template <typename TGridFunction>
2157 : number H1SemiNorm(SmartPtr<TGridFunction> spGridFct, const char* cmp, int quadOrder,
2158 : const char* subsets, ConstSmartPtr<typename H1SemiIntegrand<TGridFunction>::weight_type> weights = SPNULL)
2159 : { return H1SemiNorm(*spGridFct, cmp, quadOrder, subsets, weights); }
2160 :
2161 : template <typename TGridFunction>
2162 0 : number H1SemiNorm( SmartPtr<TGridFunction> spGridFct, const char* cmp, int quadOrder)
2163 0 : { return H1SemiNorm(*spGridFct, cmp, quadOrder, NULL); }
2164 :
2165 : template <typename TGridFunction>
2166 : number H1SemiNorm( SmartPtr<TGridFunction> spGridFct, const char* cmp, int quadOrder,
2167 : ConstSmartPtr<typename H1SemiIntegrand<TGridFunction>::weight_type> weights)
2168 : { return H1SemiNorm(*spGridFct, cmp, quadOrder, NULL, weights); }
2169 :
2170 : /// Integrand for the distance of two grid functions - evaluated in the (weighted) H1-semi norm
2171 : template <typename TGridFunction>
2172 : class H1SemiDistIntegrand : public StdIntegrand<number, TGridFunction::dim, H1SemiDistIntegrand<TGridFunction> >
2173 : {
2174 : public:
2175 : /// world dimension of grid function
2176 : static const int worldDim = TGridFunction::dim;
2177 : typedef typename H1SemiIntegrand<TGridFunction>::weight_type weight_type;
2178 :
2179 : private:
2180 : ScalarGridFunctionData<TGridFunction> m_fineData;
2181 : const int m_fineTopLevel;
2182 :
2183 : ScalarGridFunctionData<TGridFunction> m_coarseData;
2184 : const int m_coarseTopLevel;
2185 :
2186 : /// multigrid
2187 : SmartPtr<MultiGrid> m_spMG;
2188 :
2189 : /// scalar weight (optional)
2190 : ConstSmartPtr<weight_type> m_spWeight;
2191 :
2192 : public:
2193 : /// constructor
2194 : H1SemiDistIntegrand(TGridFunction& fineGridFct, size_t fineCmp,
2195 : TGridFunction& coarseGridFct, size_t coarseCmp)
2196 : : m_fineData(fineGridFct, fineCmp),
2197 : m_fineTopLevel(fineGridFct.dof_distribution()->grid_level().level()),
2198 : m_coarseData(coarseGridFct, coarseCmp),
2199 : m_coarseTopLevel(coarseGridFct.dof_distribution()->grid_level().level()),
2200 : m_spMG(fineGridFct.domain()->grid()),
2201 : m_spWeight(new ConstUserNumber<TGridFunction::dim>(1.0))
2202 : {
2203 : if(m_fineTopLevel < m_coarseTopLevel)
2204 : UG_THROW("H1SemiDiffIntegrand: fine and top level inverted.");
2205 :
2206 : if(m_fineData.domain().get() != m_coarseData.domain().get())
2207 : UG_THROW("H1SemiDiffIntegrand: grid functions defined on different domains.");
2208 : };
2209 :
2210 : /// constructor
2211 0 : H1SemiDistIntegrand(TGridFunction& fineGridFct, size_t fineCmp,
2212 : TGridFunction& coarseGridFct, size_t coarseCmp,
2213 : ConstSmartPtr<weight_type> spWeight)
2214 : : m_fineData(fineGridFct, fineCmp),
2215 0 : m_fineTopLevel(fineGridFct.dof_distribution()->grid_level().level()),
2216 : m_coarseData(coarseGridFct, coarseCmp),
2217 0 : m_coarseTopLevel(coarseGridFct.dof_distribution()->grid_level().level()),
2218 0 : m_spMG(fineGridFct.domain()->grid()),
2219 0 : m_spWeight(spWeight)
2220 : {
2221 0 : UG_COND_THROW(m_fineTopLevel < m_coarseTopLevel, "H1SemiDiffIntegrand: fine and top level inverted.");
2222 0 : UG_COND_THROW(m_fineData.domain().get() != m_coarseData.domain().get(), "H1SemiDiffIntegrand: grid functions defined on different domains.");
2223 0 : }
2224 0 : virtual ~H1SemiDistIntegrand(){}
2225 :
2226 : /// sets subset
2227 0 : virtual void set_subset(int si)
2228 : {
2229 :
2230 0 : UG_COND_THROW(!m_fineData.is_def_in_subset(si), "H1SemiDiffIntegrand: Grid function component"
2231 : <<m_fineData.fct()<<" not defined on subset "<<si);
2232 0 : UG_COND_THROW(!m_coarseData.is_def_in_subset(si), "H1SemiDiffIntegrand: Grid function component"
2233 : <<m_coarseData.fct()<<" not defined on subset "<<si);
2234 : IIntegrand<number, worldDim>::set_subset(si);
2235 0 : }
2236 :
2237 : /// \copydoc IIntegrand::values
2238 : template <int elemDim>
2239 0 : void evaluate(number vValue[],
2240 : const MathVector<worldDim> vFineGlobIP[],
2241 : GridObject* pFineElem,
2242 : const MathVector<worldDim> vCornerCoords[],
2243 : const MathVector<elemDim> vFineLocIP[],
2244 : const MathMatrix<elemDim, worldDim> vJT[],
2245 : const size_t numIP)
2246 : {
2247 : typedef typename TGridFunction::template dim_traits<elemDim>::grid_base_object Element;
2248 :
2249 : const TGridFunction &fineGridFct = m_fineData.grid_function();
2250 : const TGridFunction &coarseGridFct = m_coarseData.grid_function();
2251 :
2252 : // get coarse element
2253 : GridObject* pCoarseElem = pFineElem;
2254 0 : if(m_coarseTopLevel < m_fineTopLevel){
2255 : int parentLevel = m_spMG->get_level(pCoarseElem);
2256 0 : while(parentLevel > m_coarseTopLevel){
2257 0 : pCoarseElem = m_spMG->get_parent(pCoarseElem);
2258 : parentLevel = m_spMG->get_level(pCoarseElem);
2259 : }
2260 : }
2261 :
2262 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
2263 0 : const ReferenceObjectID fineROID = pFineElem->reference_object_id();
2264 0 : const ReferenceObjectID coarseROID = pCoarseElem->reference_object_id();
2265 :
2266 : // get corner coordinates
2267 : std::vector<MathVector<worldDim> > vCornerCoarse;
2268 0 : CollectCornerCoordinates(vCornerCoarse, *static_cast<Element*>(pCoarseElem), *m_coarseData.domain());
2269 :
2270 : // get reference Mapping
2271 : DimReferenceMapping<elemDim, worldDim>& map
2272 : = ReferenceMappingProvider::get<elemDim, worldDim>(coarseROID, vCornerCoarse);
2273 : DimReferenceMapping<elemDim, worldDim>& mapF
2274 : = ReferenceMappingProvider::get<elemDim, worldDim>(fineROID, vCornerCoords);
2275 :
2276 : std::vector<MathVector<elemDim> > vCoarseLocIP;
2277 0 : vCoarseLocIP.resize(numIP);
2278 0 : for(size_t ip = 0; ip < vCoarseLocIP.size(); ++ip) VecSet(vCoarseLocIP[ip], 0.0);
2279 0 : map.global_to_local(&vCoarseLocIP[0], vFineGlobIP, numIP);
2280 :
2281 :
2282 : // determine weights
2283 0 : std::vector<typename weight_type::data_type> elemWeights(numIP, MathMatrix<worldDim, worldDim>());
2284 : UG_ASSERT(m_spWeight.valid(), "H1SemiDistIntegrand::evaluate requires valid weights!");
2285 :
2286 : // get local solution (if required)
2287 0 : if(m_spWeight->requires_grid_fct())
2288 : {
2289 :
2290 : LocalIndices ind;
2291 : LocalVector u;
2292 : fineGridFct.indices(pFineElem, ind); // get global indices
2293 0 : u.resize(ind); // adapt local algebra
2294 0 : GetLocalVector(u, fineGridFct); // read local values of u
2295 :
2296 : // compute data
2297 : try{
2298 0 : (*m_spWeight)(&elemWeights[0], vFineGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), pFineElem,
2299 : vCornerCoords, vFineLocIP, numIP, &u, NULL);
2300 0 : } UG_CATCH_THROW("H1SemiDistIntegrand: Cannot evaluate data.");
2301 : }
2302 : else
2303 : {
2304 : // compute data
2305 : try{
2306 0 : (*m_spWeight)(&elemWeights[0], vFineGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), numIP);
2307 0 : } UG_CATCH_THROW("H1SemiDistIntegrand: Cannot evaluate data.");
2308 : }
2309 :
2310 : try{
2311 :
2312 :
2313 : // get trial space
2314 : const LocalShapeFunctionSet<elemDim>& rFineLSFS =
2315 : LocalFiniteElementProvider::get<elemDim>(fineROID, m_fineData.id());
2316 : const LocalShapeFunctionSet<elemDim>& rCoarseLSFS =
2317 : LocalFiniteElementProvider::get<elemDim>(coarseROID, m_coarseData.id());
2318 :
2319 : // get multiindices of element
2320 : std::vector<DoFIndex> vFineMI, vCoarseMI;
2321 : m_fineData.dof_indices(pFineElem, vFineMI);
2322 : m_coarseData.dof_indices(pCoarseElem, vCoarseMI);
2323 :
2324 0 : std::vector<MathVector<elemDim> > vFineLocGradient(vFineMI.size());
2325 0 : std::vector<MathVector<elemDim> > vCoarseLocGradient(vCoarseMI.size());
2326 :
2327 : // loop all integration points
2328 0 : for(size_t ip = 0; ip < numIP; ++ip)
2329 : {
2330 : // compute shape gradients at ip
2331 0 : rFineLSFS.grads(&vFineLocGradient[0], vFineLocIP[ip]);
2332 0 : rCoarseLSFS.grads(&vCoarseLocGradient[0], vCoarseLocIP[ip]);
2333 :
2334 : // compute approximated solutions at integration point
2335 : number fineSolIP = 0.0;
2336 : MathVector<elemDim> fineLocTmp(0.0);
2337 0 : for(size_t sh = 0; sh < vFineMI.size(); ++sh)
2338 : {
2339 0 : const number val = DoFRef(fineGridFct, vFineMI[sh]);
2340 0 : fineSolIP += val * rFineLSFS.shape(sh, vFineLocIP[ip]);
2341 : VecScaleAppend(fineLocTmp, val, vFineLocGradient[sh]);
2342 : }
2343 :
2344 : number coarseSolIP = 0.0;
2345 : MathVector<elemDim> coarseLocTmp(0.0);
2346 0 : for(size_t sh = 0; sh < vCoarseMI.size(); ++sh)
2347 : {
2348 0 : const number val = DoFRef(coarseGridFct, vCoarseMI[sh]);
2349 0 : coarseSolIP += val * rCoarseLSFS.shape(sh, vCoarseLocIP[ip]);
2350 : VecScaleAppend(coarseLocTmp, val, vCoarseLocGradient[sh]);
2351 : }
2352 :
2353 : // compute global gradient
2354 : MathVector<worldDim> fineGradIP;
2355 : MathMatrix<worldDim, elemDim> fineJTInv;
2356 0 : mapF.jacobian_transposed_inverse(fineJTInv, vFineLocIP[ip]);//Inverse(fineJTInv, vJT[ip]);
2357 :
2358 : #ifdef UG_DEBUG
2359 : CheckGeneralizedInverse<worldDim,elemDim>(vJT[ip], fineJTInv);
2360 : #endif //UG_DEBUG
2361 :
2362 : MatVecMult(fineGradIP, fineJTInv, fineLocTmp);
2363 :
2364 : // compute global gradient
2365 : MathVector<worldDim> coarseGradIP;
2366 : MathMatrix<worldDim, elemDim> coarseJTInv;
2367 0 : map.jacobian_transposed_inverse(coarseJTInv, vCoarseLocIP[ip]);
2368 : MatVecMult(coarseGradIP, coarseJTInv, coarseLocTmp);
2369 :
2370 : // get squared of difference
2371 : /*vValue[ip] = (coarseSolIP - fineSolIP);
2372 : vValue[ip] *= vValue[ip]; */
2373 0 : vValue[ip] = VecDistanceSq(coarseGradIP, fineGradIP, elemWeights[ip]);
2374 : }
2375 :
2376 0 : }
2377 0 : UG_CATCH_THROW("H1SemiDiffIntegrand::evaluate: trial space missing.");
2378 0 : }
2379 : };
2380 :
2381 :
2382 : //! Distance in H1 semi norm (with subset selection)
2383 : template <typename TGridFunction>
2384 : number H1SemiError2(SmartPtr<TGridFunction> spGridFct1, const char* cmp1,
2385 : SmartPtr<TGridFunction> spGridFct2, const char* cmp2,
2386 : int quadOrder, const char* subsets)
2387 : {
2388 : return GridFunctionDistance2<H1SemiDistIntegrand<TGridFunction>, TGridFunction>
2389 : (*spGridFct1, cmp1, *spGridFct2, cmp2, quadOrder, subsets);
2390 : }
2391 :
2392 : //! Distance in H1 semi norm (with subset selection)
2393 : template <typename TGridFunction>
2394 : number H1SemiError(SmartPtr<TGridFunction> spGridFct1, const char* cmp1,
2395 : SmartPtr<TGridFunction> spGridFct2, const char* cmp2,
2396 : int quadOrder, const char* subsets)
2397 : {
2398 : return sqrt(H1SemiError2(*spGridFct1, cmp1, *spGridFct2, cmp2, quadOrder, subsets));
2399 : }
2400 :
2401 : //! Distance in H1 semi norm (all subsets)
2402 : template <typename TGridFunction>
2403 : number H1SemiError(SmartPtr<TGridFunction> spGridFct1, const char* cmp1,
2404 : SmartPtr<TGridFunction> spGridFct2, const char* cmp2,
2405 : int quadOrder)
2406 : {
2407 : return H1SemiError(*spGridFct1, cmp1, *spGridFct2, cmp2, quadOrder, NULL);
2408 : }
2409 :
2410 : //! Squared distance in H1 semi norm (with select subsets & weights)
2411 : template <typename TGridFunction>
2412 0 : number H1SemiDistance2(TGridFunction& spGridFct1, const char* cmp1,
2413 : TGridFunction& spGridFct2, const char* cmp2,
2414 : int quadOrder, const char* subsets,
2415 : ConstSmartPtr<typename H1SemiDistIntegrand<TGridFunction>::weight_type> weights)
2416 : {
2417 : return GridFunctionDistance2<H1SemiDistIntegrand<TGridFunction>, TGridFunction>
2418 0 : (spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, subsets, weights);
2419 : }
2420 :
2421 : //! Distance in H1 semi norm (with select subsets & weights)
2422 : template <typename TGridFunction>
2423 : number H1SemiDistance(TGridFunction& spGridFct1, const char* cmp1,
2424 : TGridFunction& spGridFct2, const char* cmp2,
2425 : int quadOrder, const char* subsets,
2426 : ConstSmartPtr<typename H1SemiDistIntegrand<TGridFunction>::weight_type> weights)
2427 : { return sqrt(H1SemiDistance2(spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, subsets, weights)); }
2428 :
2429 : //! Squared distance in H1 semi norm (all subsets, with weights)
2430 : template <typename TGridFunction>
2431 0 : number H1SemiDistance2(TGridFunction& spGridFct1, const char* cmp1,
2432 : TGridFunction& spGridFct2, const char* cmp2,
2433 : int quadOrder, ConstSmartPtr<typename H1SemiDistIntegrand<TGridFunction>::weight_type > weights)
2434 0 : { return H1SemiDistance2(spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, NULL, weights); }
2435 :
2436 : //! Distance in H1 semi norm (all subsets, with weights)
2437 : template <typename TGridFunction>
2438 : number H1SemiDistance(TGridFunction& spGridFct1, const char* cmp1,
2439 : TGridFunction& spGridFct2, const char* cmp2,
2440 : int quadOrder, ConstSmartPtr<typename H1SemiDistIntegrand<TGridFunction>::weight_type > weights)
2441 : { return sqrt(H1SemiDistance2(spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, NULL, weights)); }
2442 :
2443 :
2444 :
2445 :
2446 : ////////////////////////////////////////////////////////////////////////////////
2447 : // H1 semi-norm Integrand
2448 : ////////////////////////////////////////////////////////////////////////////////
2449 :
2450 : //! Norm of a grid function, evaluated in (weighted) H1-semi norm
2451 : template <typename TGridFunction>
2452 : class H1EnergyIntegrand
2453 : : public StdIntegrand<number, TGridFunction::dim, H1EnergyIntegrand<TGridFunction> >
2454 : {
2455 : public:
2456 : /// world dimension of grid function
2457 : static const int worldDim = TGridFunction::dim;
2458 : typedef UserData<MathMatrix<worldDim, worldDim>, worldDim> weight_type;
2459 :
2460 : protected:
2461 : /// grid function data
2462 : ScalarGridFunctionData<TGridFunction> m_scalarData;
2463 :
2464 : /// scalar weight (optional)
2465 : ConstSmartPtr<weight_type> m_spWeight;
2466 :
2467 : public:
2468 : /// constructor
2469 0 : H1EnergyIntegrand(TGridFunction& gridFct, size_t cmp)
2470 0 : : m_scalarData(gridFct, cmp), m_spWeight(make_sp(new ConstUserMatrix<worldDim>(1.0))) {}
2471 :
2472 : /// constructor
2473 0 : H1EnergyIntegrand(TGridFunction& gridFct, size_t cmp, ConstSmartPtr<weight_type> spWeight)
2474 0 : : m_scalarData(gridFct, cmp), m_spWeight(spWeight) {}
2475 :
2476 : /// DTOR
2477 0 : virtual ~H1EnergyIntegrand() {};
2478 :
2479 : /// sets subset
2480 0 : virtual void set_subset(int si)
2481 : {
2482 :
2483 0 : UG_COND_THROW(!m_scalarData.is_def_in_subset(si), "H1EnergyIntegrand: Grid function component"
2484 : <<m_scalarData.fct()<<" not defined on subset "<<si);
2485 : IIntegrand<number, worldDim>::set_subset(si);
2486 0 : }
2487 :
2488 : /// \copydoc IIntegrand::values
2489 : template <int elemDim>
2490 0 : void evaluate(number vValue[],
2491 : const MathVector<worldDim> vGlobIP[],
2492 : GridObject* pElem,
2493 : const MathVector<worldDim> vCornerCoords[],
2494 : const MathVector<elemDim> vLocIP[],
2495 : const MathMatrix<elemDim, worldDim> vJT[],
2496 : const size_t numIP)
2497 : {
2498 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
2499 0 : const ReferenceObjectID roid = pElem->reference_object_id();
2500 : const TGridFunction &gridFct= m_scalarData.grid_function();
2501 :
2502 : DimReferenceMapping<elemDim, worldDim>& map
2503 : = ReferenceMappingProvider::get<elemDim, worldDim>(roid, vCornerCoords);
2504 :
2505 : typedef typename weight_type::data_type ipdata_type;
2506 :
2507 0 : std::vector<ipdata_type> elemWeights(numIP, MathMatrix<worldDim, worldDim>());
2508 : UG_ASSERT(m_spWeight.valid(), "H1EnergyIntegrand::evaluate requires valid weights!");
2509 :
2510 :
2511 0 : if(m_spWeight->requires_grid_fct())
2512 : {
2513 : // get local solution if needed
2514 : LocalIndices ind;
2515 : LocalVector uloc;
2516 : gridFct.indices(pElem, ind); // get global indices
2517 0 : uloc.resize(ind); // adapt local algebra
2518 0 : GetLocalVector(uloc, gridFct); // read local values of u
2519 :
2520 : // compute data
2521 : try{
2522 0 : (*m_spWeight)(&elemWeights[0], vGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), pElem,
2523 : vCornerCoords, vLocIP, numIP, &uloc, NULL);
2524 0 : } UG_CATCH_THROW("H1EnergyIntegrand: Cannot evaluate weight data.");
2525 : }
2526 : else
2527 : {
2528 : // compute data
2529 : try{
2530 0 : (*m_spWeight)(&elemWeights[0], vGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), numIP);
2531 0 : } UG_CATCH_THROW("H1EnergyIntegrand: Cannot evaluate weight data.");
2532 : }
2533 :
2534 :
2535 : try{
2536 :
2537 :
2538 :
2539 :
2540 : // get trial space
2541 : const LocalShapeFunctionSet<elemDim>& rTrialSpace =
2542 : LocalFiniteElementProvider::get<elemDim>(roid, m_scalarData.id());
2543 :
2544 : // number of dofs on element
2545 0 : const size_t num_sh = rTrialSpace.num_sh();
2546 :
2547 : // get multiindices of element
2548 : std::vector<DoFIndex> ind; // aux. index array
2549 : gridFct.dof_indices(pElem, m_scalarData.fct(), ind);
2550 :
2551 : // check multi indices
2552 0 : UG_COND_THROW(ind.size() != num_sh, "H1EnergyIntegrand::evaluate: Wrong number of multi-)indices.");
2553 :
2554 : // loop all integration points
2555 0 : std::vector<MathVector<elemDim> > vLocGradient(num_sh);
2556 0 : for(size_t ip = 0; ip < numIP; ++ip)
2557 : {
2558 : // compute shape gradients at ip
2559 0 : rTrialSpace.grads(&vLocGradient[0], vLocIP[ip]);
2560 :
2561 : // compute approximated solution at integration point
2562 : number approxSolIP = 0.0;
2563 : MathVector<elemDim> tmpVec(0.0);
2564 0 : for(size_t sh = 0; sh < num_sh; ++sh)
2565 : {
2566 : // get value at shape point (e.g. corner for P1 fct)
2567 0 : const number valSH = DoFRef(gridFct, ind[sh]);
2568 :
2569 : // add shape fct at ip * value at shape
2570 0 : approxSolIP += valSH * rTrialSpace.shape(sh, vLocIP[ip]);
2571 :
2572 : // add gradient at ip
2573 : VecScaleAppend(tmpVec, valSH, vLocGradient[sh]);
2574 : }
2575 :
2576 : // compute gradient
2577 : MathVector<worldDim> approxGradIP;
2578 : MathMatrix<worldDim, elemDim> JTInv;
2579 0 : map.jacobian_transposed_inverse(JTInv, vLocIP[ip]);//Inverse(JTInv, vJT[ip]);
2580 :
2581 : #ifdef UG_DEBUG
2582 : CheckGeneralizedInverse<worldDim,elemDim>(vJT[ip], JTInv);
2583 : #endif //UG_DEBUG
2584 :
2585 :
2586 : MatVecMult(approxGradIP, JTInv, tmpVec);
2587 :
2588 : // get norm squared
2589 : MathVector<worldDim> approxDGradIP;
2590 : MatVecMult(approxDGradIP, elemWeights[ip], approxGradIP);
2591 0 : vValue[ip] = VecTwoNormSq(approxDGradIP);
2592 : }
2593 :
2594 0 : }
2595 0 : UG_CATCH_THROW("H1EnergyIntegrand::evaluate: trial space missing.");
2596 0 : }
2597 : };
2598 :
2599 :
2600 : /// compute energy -norm of a function on the whole domain (or on selected subsets)
2601 : /**
2602 : * This function computes the integral over \f$ \| q \|^2 \f$ where the velocity is given by \f$ q:= \kappa \nabla u \f$ of a grid function u.
2603 : *
2604 : * \param[in] spGridFct grid function
2605 : * \param[in] cmp symbolic name of component function
2606 : * \param[in] time time point
2607 : * \param[in] quadOrder order of quadrature rule
2608 : * \param[in] subsets subsets, where to compute (OPTIONAL)
2609 : * \param[in] weights element-wise matrix-valued weights kappa (OPTIONAL)
2610 : * \returns number l2-norm
2611 : */
2612 : template <typename TGridFunction>
2613 0 : number H1EnergyNorm2(TGridFunction& gridFct, const char* cmp, int quadOrder, const char* subsets=NULL,
2614 : ConstSmartPtr<typename H1SemiIntegrand<TGridFunction>::weight_type> weights = SPNULL)
2615 : {
2616 : // get function id of name
2617 : const size_t fct = gridFct.fct_id_by_name(cmp);
2618 :
2619 : // check that function exists
2620 0 : UG_COND_THROW(fct >= gridFct.num_fct(), "H1SemiNorm: Function space does not contain"
2621 : " a function with name " << cmp << ".");
2622 0 : if (weights.invalid()) {
2623 0 : H1EnergyIntegrand<TGridFunction> integrand(gridFct, fct);
2624 0 : return IntegrateSubsets(integrand, gridFct, subsets, quadOrder);
2625 : } else {
2626 0 : H1EnergyIntegrand<TGridFunction> integrand(gridFct, fct, weights);
2627 0 : return IntegrateSubsets(integrand, gridFct, subsets, quadOrder);
2628 : }
2629 : }
2630 : template <typename TGridFunction>
2631 : number H1EnergyNorm(TGridFunction& gridFct, const char* cmp, int quadOrder, const char* subsets=NULL,
2632 : ConstSmartPtr<typename H1SemiIntegrand<TGridFunction>::weight_type> weights = SPNULL)
2633 : {
2634 : return (sqrt(H1EnergyNorm2(gridFct, cmp, quadOrder, subsets, weights)));
2635 : }
2636 :
2637 : template <typename TGridFunction>
2638 : number H1EnergyNorm(SmartPtr<TGridFunction> spGridFct, const char* cmp, int quadOrder,
2639 : const char* subsets, ConstSmartPtr<typename H1SemiIntegrand<TGridFunction>::weight_type> weights = SPNULL)
2640 : { return H1EnergyNorm(*spGridFct, cmp, quadOrder, subsets, weights); }
2641 :
2642 : template <typename TGridFunction>
2643 : number H1EnergyNorm( SmartPtr<TGridFunction> spGridFct, const char* cmp, int quadOrder)
2644 : { return H1EnergyNorm(spGridFct, cmp, quadOrder, NULL); }
2645 :
2646 : template <typename TGridFunction>
2647 : number H1EnergyNorm( SmartPtr<TGridFunction> spGridFct, const char* cmp, int quadOrder,
2648 : ConstSmartPtr<typename H1SemiIntegrand<TGridFunction>::weight_type> weights)
2649 : { return H1EnergyNorm(spGridFct, cmp, quadOrder, NULL, weights); }
2650 :
2651 :
2652 :
2653 :
2654 : /// Integrand for the distance of two grid functions - evaluated in the norm \f$ |D \nabla u|^2 \f$
2655 : template <typename TGridFunction>
2656 : class H1EnergyDistIntegrand
2657 : : public StdIntegrand<number, TGridFunction::dim, H1EnergyDistIntegrand<TGridFunction> >
2658 : {
2659 : public:
2660 : /// world dimension of grid function
2661 : static const int worldDim = TGridFunction::dim;
2662 : typedef typename H1SemiIntegrand<TGridFunction>::weight_type weight_type;
2663 :
2664 : private:
2665 : ScalarGridFunctionData<TGridFunction> m_fineData;
2666 : const int m_fineTopLevel;
2667 :
2668 : ScalarGridFunctionData<TGridFunction> m_coarseData;
2669 : const int m_coarseTopLevel;
2670 :
2671 : /// multigrid
2672 : SmartPtr<MultiGrid> m_spMG;
2673 :
2674 : /// scalar weight (optional)
2675 : ConstSmartPtr<weight_type> m_spWeight;
2676 :
2677 : public:
2678 : /// constructor
2679 : H1EnergyDistIntegrand(TGridFunction& fineGridFct, size_t fineCmp,
2680 : TGridFunction& coarseGridFct, size_t coarseCmp)
2681 : : m_fineData(fineGridFct, fineCmp),
2682 : m_fineTopLevel(fineGridFct.dof_distribution()->grid_level().level()),
2683 : m_coarseData(coarseGridFct, coarseCmp),
2684 : m_coarseTopLevel(coarseGridFct.dof_distribution()->grid_level().level()),
2685 : m_spMG(fineGridFct.domain()->grid()),
2686 : m_spWeight(new ConstUserNumber<TGridFunction::dim>(1.0))
2687 : {
2688 : if(m_fineTopLevel < m_coarseTopLevel)
2689 : UG_THROW("H1EnergyDistIntegrand: fine and top level inverted.");
2690 :
2691 : if(m_fineData.domain().get() != m_coarseData.domain().get())
2692 : UG_THROW("H1EnergyDistIntegrand: grid functions defined on different domains.");
2693 : };
2694 :
2695 : /// constructor
2696 0 : H1EnergyDistIntegrand(TGridFunction& fineGridFct, size_t fineCmp,
2697 : TGridFunction& coarseGridFct, size_t coarseCmp,
2698 : ConstSmartPtr<weight_type> spWeight)
2699 : : m_fineData(fineGridFct, fineCmp),
2700 0 : m_fineTopLevel(fineGridFct.dof_distribution()->grid_level().level()),
2701 : m_coarseData(coarseGridFct, coarseCmp),
2702 0 : m_coarseTopLevel(coarseGridFct.dof_distribution()->grid_level().level()),
2703 0 : m_spMG(fineGridFct.domain()->grid()),
2704 0 : m_spWeight(spWeight)
2705 : {
2706 0 : if(m_fineTopLevel < m_coarseTopLevel)
2707 0 : UG_THROW("H1EnergyDistIntegrand: fine and top level inverted.");
2708 :
2709 0 : if(m_fineData.domain().get() != m_coarseData.domain().get())
2710 0 : UG_THROW("H1EnergyDistIntegrand: grid functions defined on different domains.");
2711 0 : }
2712 0 : virtual ~H1EnergyDistIntegrand(){}
2713 :
2714 : /// sets subset
2715 0 : virtual void set_subset(int si)
2716 : {
2717 0 : if(!m_fineData.is_def_in_subset(si))
2718 0 : UG_THROW("H1EnergyDistIntegrand: Grid function component"
2719 : <<m_fineData.fct()<<" not defined on subset "<<si);
2720 0 : if(!m_coarseData.is_def_in_subset(si))
2721 0 : UG_THROW("H1EnergyDistIntegrand: Grid function component"
2722 : <<m_coarseData.fct()<<" not defined on subset "<<si);
2723 : IIntegrand<number, worldDim>::set_subset(si);
2724 0 : }
2725 :
2726 : /// \copydoc IIntegrand::values
2727 : template <int elemDim>
2728 0 : void evaluate(number vValue[],
2729 : const MathVector<worldDim> vFineGlobIP[],
2730 : GridObject* pFineElem,
2731 : const MathVector<worldDim> vCornerCoords[],
2732 : const MathVector<elemDim> vFineLocIP[],
2733 : const MathMatrix<elemDim, worldDim> vJT[],
2734 : const size_t numIP)
2735 : {
2736 : typedef typename TGridFunction::template dim_traits<elemDim>::grid_base_object Element;
2737 :
2738 : const TGridFunction &fineGridFct = m_fineData.grid_function();
2739 : const TGridFunction &coarseGridFct = m_coarseData.grid_function();
2740 :
2741 : // get coarse element
2742 : GridObject* pCoarseElem = pFineElem;
2743 0 : if(m_coarseTopLevel < m_fineTopLevel){
2744 : int parentLevel = m_spMG->get_level(pCoarseElem);
2745 0 : while(parentLevel > m_coarseTopLevel){
2746 0 : pCoarseElem = m_spMG->get_parent(pCoarseElem);
2747 : parentLevel = m_spMG->get_level(pCoarseElem);
2748 : }
2749 : }
2750 :
2751 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
2752 0 : const ReferenceObjectID fineROID = pFineElem->reference_object_id();
2753 0 : const ReferenceObjectID coarseROID = pCoarseElem->reference_object_id();
2754 :
2755 : // get corner coordinates
2756 : std::vector<MathVector<worldDim> > vCornerCoarse;
2757 0 : CollectCornerCoordinates(vCornerCoarse, *static_cast<Element*>(pCoarseElem), *m_coarseData.domain());
2758 :
2759 : // get reference Mapping
2760 : DimReferenceMapping<elemDim, worldDim>& map
2761 : = ReferenceMappingProvider::get<elemDim, worldDim>(coarseROID, vCornerCoarse);
2762 : DimReferenceMapping<elemDim, worldDim>& mapF
2763 : = ReferenceMappingProvider::get<elemDim, worldDim>(fineROID, vCornerCoords);
2764 :
2765 : std::vector<MathVector<elemDim> > vCoarseLocIP;
2766 0 : vCoarseLocIP.resize(numIP);
2767 0 : for(size_t ip = 0; ip < vCoarseLocIP.size(); ++ip) VecSet(vCoarseLocIP[ip], 0.0);
2768 0 : map.global_to_local(&vCoarseLocIP[0], vFineGlobIP, numIP);
2769 :
2770 :
2771 : // determine weights
2772 0 : std::vector<typename weight_type::data_type> elemWeights(numIP, MathMatrix<worldDim, worldDim>());
2773 : UG_ASSERT(m_spWeight.valid(), "H1SemiDistIntegrand::evaluate requires valid weights!");
2774 :
2775 : // get local solution if needed
2776 0 : if(m_spWeight->requires_grid_fct())
2777 : {
2778 :
2779 : LocalIndices ind;
2780 : LocalVector u;
2781 : fineGridFct.indices(pFineElem, ind); // get global indices
2782 0 : u.resize(ind); // adapt local algebra
2783 0 : GetLocalVector(u, fineGridFct); // read local values of u
2784 :
2785 : // compute data
2786 : try{
2787 0 : (*m_spWeight)(&elemWeights[0], vFineGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), pFineElem,
2788 : vCornerCoords, vFineLocIP, numIP, &u, NULL);
2789 0 : } UG_CATCH_THROW("H1SemiDistIntegrand: Cannot evaluate data.");
2790 : }
2791 : else
2792 : {
2793 : // compute data
2794 : try{
2795 0 : (*m_spWeight)(&elemWeights[0], vFineGlobIP, 0.0, IIntegrand<number, worldDim>::subset(), numIP);
2796 0 : } UG_CATCH_THROW("H1SemiDistIntegrand: Cannot evaluate data.");
2797 : }
2798 :
2799 : try{
2800 :
2801 :
2802 : // get trial space
2803 : const LocalShapeFunctionSet<elemDim>& rFineLSFS =
2804 : LocalFiniteElementProvider::get<elemDim>(fineROID, m_fineData.id());
2805 : const LocalShapeFunctionSet<elemDim>& rCoarseLSFS =
2806 : LocalFiniteElementProvider::get<elemDim>(coarseROID, m_coarseData.id());
2807 :
2808 : // get multiindices of element
2809 : std::vector<DoFIndex> vFineMI, vCoarseMI;
2810 : m_fineData.dof_indices(pFineElem, vFineMI);
2811 : m_coarseData.dof_indices(pCoarseElem, vCoarseMI);
2812 :
2813 0 : std::vector<MathVector<elemDim> > vFineLocGradient(vFineMI.size());
2814 0 : std::vector<MathVector<elemDim> > vCoarseLocGradient(vCoarseMI.size());
2815 :
2816 : // loop all integration points
2817 0 : for(size_t ip = 0; ip < numIP; ++ip)
2818 : {
2819 : // compute shape gradients at ip
2820 0 : rFineLSFS.grads(&vFineLocGradient[0], vFineLocIP[ip]);
2821 0 : rCoarseLSFS.grads(&vCoarseLocGradient[0], vCoarseLocIP[ip]);
2822 :
2823 : // compute approximated solutions at integration point
2824 : number fineSolIP = 0.0;
2825 : MathVector<elemDim> fineLocTmp(0.0);
2826 0 : for(size_t sh = 0; sh < vFineMI.size(); ++sh)
2827 : {
2828 0 : const number val = DoFRef(fineGridFct, vFineMI[sh]);
2829 0 : fineSolIP += val * rFineLSFS.shape(sh, vFineLocIP[ip]);
2830 : VecScaleAppend(fineLocTmp, val, vFineLocGradient[sh]);
2831 : }
2832 :
2833 : number coarseSolIP = 0.0;
2834 : MathVector<elemDim> coarseLocTmp(0.0);
2835 0 : for(size_t sh = 0; sh < vCoarseMI.size(); ++sh)
2836 : {
2837 0 : const number val = DoFRef(coarseGridFct, vCoarseMI[sh]);
2838 0 : coarseSolIP += val * rCoarseLSFS.shape(sh, vCoarseLocIP[ip]);
2839 : VecScaleAppend(coarseLocTmp, val, vCoarseLocGradient[sh]);
2840 : }
2841 :
2842 : // compute global D*gradient
2843 : MathVector<worldDim> fineGradIP;
2844 : MathMatrix<worldDim, elemDim> fineJTInv;
2845 0 : mapF.jacobian_transposed_inverse(fineJTInv, vFineLocIP[ip]);//Inverse(fineJTInv, vJT[ip]);
2846 :
2847 : #ifdef UG_DEBUG
2848 : CheckGeneralizedInverse<worldDim,elemDim>(vJT[ip], fineJTInv);
2849 : #endif //UG_DEBUG
2850 :
2851 :
2852 : MatVecMult(fineGradIP, fineJTInv, fineLocTmp);
2853 : MathVector<worldDim> fineWorldLocTmp(0.0);
2854 : MatVecMult(fineWorldLocTmp, elemWeights[ip], fineGradIP);
2855 :
2856 : // compute global D*gradient
2857 : MathVector<worldDim> coarseGradIP;
2858 : MathMatrix<worldDim, elemDim> coarseJTInv;
2859 0 : map.jacobian_transposed_inverse(coarseJTInv, vCoarseLocIP[ip]);
2860 : MatVecMult(coarseGradIP, coarseJTInv, coarseLocTmp);
2861 : MathVector<worldDim> coarseWorldLocTmp(0.0);
2862 : MatVecMult(coarseWorldLocTmp, elemWeights[ip], coarseGradIP);
2863 :
2864 : // get squared difference
2865 0 : vValue[ip] = VecDistanceSq(fineWorldLocTmp, coarseWorldLocTmp);
2866 : //vValue[ip] = VecDistanceSq(fineGradIP, coarseGradIP, elemWeights[ip]);
2867 : }
2868 :
2869 0 : }
2870 0 : UG_CATCH_THROW("H1EnergyDiffIntegrand::evaluate: trial space missing.");
2871 0 : }
2872 : };
2873 :
2874 :
2875 : //! Squared distance in H1 semi norm (with select subsets & weights)
2876 : template <typename TGridFunction>
2877 0 : number H1EnergyDistance2(TGridFunction& spGridFct1, const char* cmp1,
2878 : TGridFunction& spGridFct2, const char* cmp2,
2879 : int quadOrder, const char* subsets,
2880 : ConstSmartPtr<typename H1SemiDistIntegrand<TGridFunction>::weight_type> weights)
2881 : {
2882 : return GridFunctionDistance2<H1EnergyDistIntegrand<TGridFunction>, TGridFunction>
2883 0 : (spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, subsets, weights);
2884 : }
2885 :
2886 : //! Distance in H1 semi norm (with select subsets & weights)
2887 : template <typename TGridFunction>
2888 : number H1EnergyDistance(TGridFunction& spGridFct1, const char* cmp1,
2889 : TGridFunction& spGridFct2, const char* cmp2,
2890 : int quadOrder, const char* subsets,
2891 : ConstSmartPtr<typename H1SemiDistIntegrand<TGridFunction>::weight_type> weights)
2892 : { return sqrt(H1EnergyDistance2(spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, subsets, weights)); }
2893 :
2894 : //! Squared distance in H1 semi norm (all subsets, with weights)
2895 : template <typename TGridFunction>
2896 : number H1EnergyDistance2(TGridFunction& spGridFct1, const char* cmp1,
2897 : TGridFunction& spGridFct2, const char* cmp2,
2898 : int quadOrder, ConstSmartPtr<typename H1SemiDistIntegrand<TGridFunction>::weight_type > weights)
2899 : { return H1EnergyDistance2(spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, NULL, weights); }
2900 :
2901 : //! Distance in H1 semi norm (all subsets, with weights)
2902 : template <typename TGridFunction>
2903 : number H1EnergyDistance(TGridFunction& spGridFct1, const char* cmp1,
2904 : TGridFunction& spGridFct2, const char* cmp2,
2905 : int quadOrder, ConstSmartPtr<typename H1SemiDistIntegrand<TGridFunction>::weight_type > weights)
2906 : { return sqrt(H1EnergyDistance2(spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, NULL, weights)); }
2907 :
2908 :
2909 :
2910 :
2911 : ////////////////////////////////////////////////////////////////////////////////
2912 : // H1 norm integrand
2913 : ////////////////////////////////////////////////////////////////////////////////
2914 : template <typename TGridFunction>
2915 : class H1NormIntegrand
2916 : : public StdIntegrand<number, TGridFunction::dim, H1NormIntegrand<TGridFunction> >
2917 : {
2918 : public:
2919 : /// world dimension of grid function
2920 : static const int worldDim = TGridFunction::dim;
2921 :
2922 : private:
2923 : ScalarGridFunctionData<TGridFunction> m_scalarData;
2924 :
2925 : public:
2926 : /// CTOR
2927 0 : H1NormIntegrand(TGridFunction& gridFct, size_t cmp)
2928 0 : : m_scalarData(gridFct, cmp) {}
2929 :
2930 : /// DTOR
2931 0 : virtual ~H1NormIntegrand() {}
2932 :
2933 : /// sets subset
2934 0 : virtual void set_subset(int si)
2935 : {
2936 0 : if(!m_scalarData.is_def_in_subset(si))
2937 0 : UG_THROW("H1Norm: Grid function component"
2938 : <<m_scalarData.fct()<<" not defined on subset "<<si);
2939 : IIntegrand<number, worldDim>::set_subset(si);
2940 0 : }
2941 :
2942 : /// \copydoc IIntegrand::values
2943 : template <int elemDim>
2944 0 : void evaluate(number vValue[],
2945 : const MathVector<worldDim> vGlobIP[],
2946 : GridObject* pElem,
2947 : const MathVector<worldDim> vCornerCoords[],
2948 : const MathVector<elemDim> vLocIP[],
2949 : const MathMatrix<elemDim, worldDim> vJT[],
2950 : const size_t numIP)
2951 : {
2952 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
2953 0 : const ReferenceObjectID roid = pElem->reference_object_id();
2954 :
2955 : DimReferenceMapping<elemDim, worldDim>& map
2956 : = ReferenceMappingProvider::get<elemDim, worldDim>(roid, vCornerCoords);
2957 :
2958 : try{
2959 : // get trial space
2960 : const LocalShapeFunctionSet<elemDim>& rTrialSpace =
2961 : LocalFiniteElementProvider::get<elemDim>(roid, m_scalarData.id());
2962 :
2963 : // number of dofs on element
2964 0 : const size_t num_sh = rTrialSpace.num_sh();
2965 :
2966 : // get multiindices of element
2967 : std::vector<DoFIndex> ind; // aux. index array
2968 : m_scalarData.dof_indices(pElem, ind);
2969 :
2970 : // check multi indices
2971 0 : if(ind.size() != num_sh)
2972 0 : UG_THROW("H1ErrorIntegrand::evaluate: Wrong number of"
2973 : " multi indices.");
2974 :
2975 : // loop all integration points
2976 0 : std::vector<MathVector<elemDim> > vLocGradient(num_sh);
2977 0 : for(size_t ip = 0; ip < numIP; ++ip)
2978 : {
2979 :
2980 : // compute shape gradients at ip
2981 0 : rTrialSpace.grads(&vLocGradient[0], vLocIP[ip]);
2982 :
2983 : // compute approximated solution at integration point
2984 : number approxSolIP = 0.0;
2985 : MathVector<elemDim> locTmp; VecSet(locTmp, 0.0);
2986 0 : for(size_t sh = 0; sh < num_sh; ++sh)
2987 : {
2988 : // get value at shape point (e.g. corner for P1 fct)
2989 0 : const number valSH = DoFRef(m_scalarData.grid_function(), ind[sh]);
2990 :
2991 : // add shape fct at ip * value at shape
2992 0 : approxSolIP += valSH * rTrialSpace.shape(sh, vLocIP[ip]);
2993 :
2994 : // add gradient at ip
2995 : VecScaleAppend(locTmp, valSH, vLocGradient[sh]);
2996 : }
2997 :
2998 : // compute global gradient
2999 : MathVector<worldDim> approxGradIP;
3000 : MathMatrix<worldDim, elemDim> JTInv;
3001 0 : map.jacobian_transposed_inverse(JTInv, vLocIP[ip]);//RightInverse(JTInv, vJT[ip]);
3002 :
3003 : #ifdef UG_DEBUG
3004 : CheckGeneralizedInverse<worldDim,elemDim>(vJT[ip], JTInv);
3005 : #endif //UG_DEBUG
3006 :
3007 : MatVecMult(approxGradIP, JTInv, locTmp);
3008 :
3009 : // get squared of difference
3010 0 : vValue[ip] = approxSolIP * approxSolIP;
3011 0 : vValue[ip] += VecDot(approxGradIP, approxGradIP);
3012 : }
3013 :
3014 0 : }
3015 0 : UG_CATCH_THROW("H1SemiNormFuncIntegrand::evaluate: trial space missing.");
3016 0 : }
3017 : };
3018 :
3019 :
3020 :
3021 :
3022 : template <typename TGridFunction>
3023 0 : number H1Norm2(TGridFunction& u, const char* cmp,
3024 : int quadOrder, const char* subsets=NULL)
3025 : {
3026 : // get function id of name
3027 : const size_t fct = u.fct_id_by_name(cmp);
3028 :
3029 : // check that function exists
3030 0 : if(fct >= u.num_fct())
3031 0 : UG_THROW("H1Norm: Function space does not contain"
3032 : " a function with name " << cmp << ".");
3033 :
3034 : H1NormIntegrand<TGridFunction> spIntegrand(u, fct);
3035 0 : return IntegrateSubsets(spIntegrand, u, subsets, quadOrder);
3036 : }
3037 :
3038 :
3039 : template <typename TGridFunction>
3040 : number H1Norm(TGridFunction& u, const char* cmp,
3041 : int quadOrder, const char* subsets=NULL)
3042 : {
3043 0 : return sqrt(H1Norm2(u, cmp, quadOrder,subsets));
3044 : }
3045 :
3046 : template <typename TGridFunction>
3047 0 : number H1Norm(SmartPtr<TGridFunction> spGridFct, const char* cmp,
3048 : int quadOrder, const char* subsets)
3049 : {
3050 0 : return H1Norm(*spGridFct, cmp, quadOrder, subsets);
3051 : }
3052 :
3053 : template <typename TGridFunction>
3054 0 : number H1Norm(SmartPtr<TGridFunction> spGridFct, const char* cmp, int quadOrder)
3055 : {
3056 0 : return H1Norm(*spGridFct, cmp, quadOrder, NULL);
3057 : }
3058 :
3059 :
3060 :
3061 : /// Integrand for the distance of two grid functions - evaluated in the H1 norm
3062 : template <typename TGridFunction>
3063 : class H1DistIntegrand
3064 : : public StdIntegrand<number, TGridFunction::dim, H1DistIntegrand<TGridFunction> >
3065 : {
3066 : public:
3067 : /// world dimension of grid function
3068 : static const int worldDim = TGridFunction::dim;
3069 :
3070 : private:
3071 : ScalarGridFunctionData<TGridFunction> m_fineData;
3072 : const int m_fineTopLevel;
3073 :
3074 : ScalarGridFunctionData<TGridFunction> m_coarseData;
3075 : const int m_coarseTopLevel;
3076 :
3077 : /// multigrid
3078 : SmartPtr<MultiGrid> m_spMG;
3079 :
3080 : public:
3081 : /// constructor (1 is fine grid function)
3082 0 : H1DistIntegrand(TGridFunction& fineGridFct, size_t fineCmp,
3083 : TGridFunction& coarseGridFct, size_t coarseCmp)
3084 : : m_fineData(fineGridFct, fineCmp),
3085 0 : m_fineTopLevel(fineGridFct.dof_distribution()->grid_level().level()),
3086 : m_coarseData(coarseGridFct, coarseCmp),
3087 0 : m_coarseTopLevel(coarseGridFct.dof_distribution()->grid_level().level()),
3088 0 : m_spMG(fineGridFct.domain()->grid())
3089 : {
3090 0 : if(m_fineTopLevel < m_coarseTopLevel)
3091 0 : UG_THROW("H1DiffIntegrand: fine and top level inverted.");
3092 :
3093 0 : if(fineGridFct.domain().get() !=
3094 0 : coarseGridFct.domain().get())
3095 0 : UG_THROW("H1DiffIntegrand: grid functions defined on different domains.");
3096 0 : };
3097 :
3098 : /// DTOR
3099 0 : virtual ~H1DistIntegrand() {};
3100 :
3101 : /// sets subset
3102 0 : virtual void set_subset(int si)
3103 : {
3104 0 : if(!m_fineData.is_def_in_subset(si))
3105 0 : UG_THROW("H1DiffIntegrand: Grid function component"
3106 : <<m_fineData.fct()<<" not defined on subset "<<si);
3107 0 : if(!m_coarseData.is_def_in_subset(si))
3108 0 : UG_THROW("H1DiffIntegrand: Grid function component"
3109 : <<m_coarseData.fct()<<" not defined on subset "<<si);
3110 : IIntegrand<number, worldDim>::set_subset(si);
3111 0 : }
3112 :
3113 : /// \copydoc IIntegrand::values
3114 : template <int elemDim>
3115 0 : void evaluate(number vValue[],
3116 : const MathVector<worldDim> vFineGlobIP[],
3117 : GridObject* pFineElem,
3118 : const MathVector<worldDim> vCornerCoords[],
3119 : const MathVector<elemDim> vFineLocIP[],
3120 : const MathMatrix<elemDim, worldDim> vJT[],
3121 : const size_t numIP)
3122 : {
3123 : typedef typename TGridFunction::template dim_traits<elemDim>::grid_base_object Element;
3124 :
3125 : // get coarse element
3126 : GridObject* pCoarseElem = pFineElem;
3127 0 : if(m_coarseTopLevel < m_fineTopLevel){
3128 : int parentLevel = m_spMG->get_level(pCoarseElem);
3129 0 : while(parentLevel > m_coarseTopLevel){
3130 0 : pCoarseElem = m_spMG->get_parent(pCoarseElem);
3131 : parentLevel = m_spMG->get_level(pCoarseElem);
3132 : }
3133 : }
3134 :
3135 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
3136 0 : const ReferenceObjectID fineROID = pFineElem->reference_object_id();
3137 0 : const ReferenceObjectID coarseROID = pCoarseElem->reference_object_id();
3138 :
3139 : // get corner coordinates
3140 : std::vector<MathVector<worldDim> > vCornerCoarse;
3141 0 : CollectCornerCoordinates(vCornerCoarse, *static_cast<Element*>(pCoarseElem), *m_coarseData.domain());
3142 :
3143 : // get Reference Mapping
3144 : DimReferenceMapping<elemDim, worldDim>& map
3145 : = ReferenceMappingProvider::get<elemDim, worldDim>(coarseROID, vCornerCoarse);
3146 : DimReferenceMapping<elemDim, worldDim>& mapF
3147 : = ReferenceMappingProvider::get<elemDim, worldDim>(fineROID, vCornerCoords);
3148 :
3149 : std::vector<MathVector<elemDim> > vCoarseLocIP;
3150 0 : vCoarseLocIP.resize(numIP);
3151 0 : for(size_t ip = 0; ip < vCoarseLocIP.size(); ++ip) VecSet(vCoarseLocIP[ip], 0.0);
3152 0 : map.global_to_local(&vCoarseLocIP[0], vFineGlobIP, numIP);
3153 :
3154 : try{
3155 : // get trial space
3156 : const LocalShapeFunctionSet<elemDim>& rFineLSFS =
3157 : LocalFiniteElementProvider::get<elemDim>(fineROID, m_fineData.id());
3158 : const LocalShapeFunctionSet<elemDim>& rCoarseLSFS =
3159 : LocalFiniteElementProvider::get<elemDim>(coarseROID, m_coarseData.id());
3160 :
3161 : // get multiindices of element
3162 : std::vector<DoFIndex> vFineMI, vCoarseMI;
3163 : m_fineData.dof_indices(pFineElem, vFineMI);
3164 : m_coarseData.dof_indices(pCoarseElem, vCoarseMI);
3165 :
3166 0 : std::vector<MathVector<elemDim> > vFineLocGradient(vFineMI.size());
3167 0 : std::vector<MathVector<elemDim> > vCoarseLocGradient(vCoarseMI.size());
3168 :
3169 : // loop all integration points
3170 0 : for(size_t ip = 0; ip < numIP; ++ip)
3171 : {
3172 : // compute shape gradients at ip
3173 0 : rFineLSFS.grads(&vFineLocGradient[0], vFineLocIP[ip]);
3174 0 : rCoarseLSFS.grads(&vCoarseLocGradient[0], vCoarseLocIP[ip]);
3175 :
3176 : // compute approximated solution at integration point
3177 : number fineSolIP = 0.0;
3178 : MathVector<elemDim> fineLocTmp; VecSet(fineLocTmp, 0.0);
3179 0 : for(size_t sh = 0; sh < vFineMI.size(); ++sh)
3180 : {
3181 0 : const number val = DoFRef(m_fineData.grid_function(), vFineMI[sh]);
3182 0 : fineSolIP += val * rFineLSFS.shape(sh, vFineLocIP[ip]);
3183 : VecScaleAppend(fineLocTmp, val, vFineLocGradient[sh]);
3184 : }
3185 : number coarseSolIP = 0.0;
3186 : MathVector<elemDim> coarseLocTmp; VecSet(coarseLocTmp, 0.0);
3187 0 : for(size_t sh = 0; sh < vCoarseMI.size(); ++sh)
3188 : {
3189 0 : const number val = DoFRef(m_coarseData.grid_function(), vCoarseMI[sh]);
3190 0 : coarseSolIP += val * rCoarseLSFS.shape(sh, vCoarseLocIP[ip]);
3191 : VecScaleAppend(coarseLocTmp, val, vCoarseLocGradient[sh]);
3192 : }
3193 :
3194 : // compute global gradient
3195 : MathVector<worldDim> fineGradIP;
3196 : MathMatrix<worldDim, elemDim> fineJTInv;
3197 0 : mapF.jacobian_transposed_inverse(fineJTInv, vFineLocIP[ip]);//RightInverse(fineJTInv, vJT[ip]);
3198 :
3199 : #ifdef UG_DEBUG
3200 : CheckGeneralizedInverse<worldDim,elemDim>(vJT[ip], fineJTInv);
3201 : #endif //UG_DEBUG
3202 :
3203 :
3204 : MatVecMult(fineGradIP, fineJTInv, fineLocTmp);
3205 :
3206 : // compute global gradient
3207 : MathVector<worldDim> coarseGradIP;
3208 : MathMatrix<worldDim, elemDim> coarseJTInv;
3209 0 : map.jacobian_transposed_inverse(coarseJTInv, vCoarseLocIP[ip]);
3210 : MatVecMult(coarseGradIP, coarseJTInv, coarseLocTmp);
3211 :
3212 : // get squared of difference
3213 0 : vValue[ip] = (coarseSolIP - fineSolIP);
3214 0 : vValue[ip] *= vValue[ip];
3215 0 : vValue[ip] += VecDistanceSq(coarseGradIP, fineGradIP);
3216 : }
3217 :
3218 0 : }
3219 0 : UG_CATCH_THROW("H1DiffIntegrand::evaluate: trial space missing.");
3220 0 : }
3221 : };
3222 :
3223 :
3224 :
3225 : template <typename TGridFunction>
3226 : number H1Distance2(TGridFunction& spGridFct1, const char* cmp1,
3227 : TGridFunction& spGridFct2, const char* cmp2,
3228 : int quadOrder, const char* subsets=NULL)
3229 : {
3230 : return GridFunctionDistance2<H1DistIntegrand<TGridFunction>, TGridFunction>
3231 0 : (spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, subsets);
3232 : }
3233 :
3234 :
3235 : template <typename TGridFunction>
3236 : number H1Distance(TGridFunction& spGridFct1, const char* cmp1,
3237 : TGridFunction& spGridFct2, const char* cmp2,
3238 : int quadOrder, const char* subsets=NULL)
3239 : {
3240 0 : return sqrt(H1Distance2(spGridFct1, cmp1, spGridFct2, cmp2, quadOrder, subsets));
3241 : }
3242 : /// for lua shell
3243 : template <typename TGridFunction>
3244 0 : number H1Error(SmartPtr<TGridFunction> spGridFct1, const char* cmp1,
3245 : SmartPtr<TGridFunction> spGridFct2, const char* cmp2,
3246 : int quadOrder, const char* subsets)
3247 : {
3248 0 : return H1Distance(*spGridFct1, cmp1, *spGridFct2, cmp2, quadOrder, subsets);
3249 : }
3250 :
3251 : /// for lua shell
3252 : template <typename TGridFunction>
3253 0 : number H1Error(SmartPtr<TGridFunction> spGridFct1, const char* cmp1,
3254 : SmartPtr<TGridFunction> spGridFct2, const char* cmp2,
3255 : int quadOrder)
3256 : {
3257 0 : return H1Distance(*spGridFct1, cmp1, *spGridFct2, cmp2, quadOrder, NULL);
3258 : }
3259 :
3260 :
3261 : ////////////////////////////////////////////////////////////////////////////////
3262 : // Standard Integrand
3263 : ////////////////////////////////////////////////////////////////////////////////
3264 :
3265 : template <typename TGridFunction>
3266 : class StdFuncIntegrand
3267 : : public StdIntegrand<number, TGridFunction::dim, StdFuncIntegrand<TGridFunction> >
3268 : {
3269 : public:
3270 : // world dimension of grid function
3271 : static const int worldDim = TGridFunction::dim;
3272 :
3273 : private:
3274 : // grid function
3275 : TGridFunction* m_pGridFct;
3276 :
3277 : // component of function
3278 : const size_t m_fct;
3279 :
3280 : public:
3281 : /// constructor
3282 0 : StdFuncIntegrand(TGridFunction* pGridFct, size_t cmp)
3283 0 : : m_pGridFct(pGridFct), m_fct(cmp)
3284 : {};
3285 :
3286 0 : virtual ~StdFuncIntegrand(){}
3287 :
3288 : /// sets subset
3289 0 : virtual void set_subset(int si)
3290 : {
3291 0 : if(!m_pGridFct->is_def_in_subset(m_fct, si))
3292 0 : UG_THROW("L2ErrorIntegrand: Grid function component"
3293 : <<m_fct<<" not defined on subset "<<si);
3294 : IIntegrand<number, worldDim>::set_subset(si);
3295 0 : }
3296 :
3297 : /// \copydoc IIntegrand::values
3298 : template <int elemDim>
3299 0 : void evaluate(number vValue[],
3300 : const MathVector<worldDim> vGlobIP[],
3301 : GridObject* pElem,
3302 : const MathVector<worldDim> vCornerCoords[],
3303 : const MathVector<elemDim> vLocIP[],
3304 : const MathMatrix<elemDim, worldDim> vJT[],
3305 : const size_t numIP)
3306 : {
3307 : // get reference object id (i.e. Triangle, Quadrilateral, Tetrahedron, ...)
3308 0 : ReferenceObjectID roid = (ReferenceObjectID) pElem->reference_object_id();
3309 :
3310 0 : const LFEID m_id = m_pGridFct->local_finite_element_id(m_fct);
3311 :
3312 : try{
3313 : // get trial space
3314 : const LocalShapeFunctionSet<elemDim>& rTrialSpace =
3315 : LocalFiniteElementProvider::get<elemDim>(roid, m_id);
3316 :
3317 : // number of dofs on element
3318 0 : const size_t num_sh = rTrialSpace.num_sh();
3319 :
3320 : // get multiindices of element
3321 :
3322 : std::vector<DoFIndex> ind; // aux. index array
3323 0 : m_pGridFct->dof_indices(pElem, m_fct, ind);
3324 :
3325 : // check multi indices
3326 0 : UG_COND_THROW(ind.size() != num_sh,
3327 : "StdFuncIntegrand::evaluate: Wrong number of multi indices.");
3328 :
3329 : // loop all integration points
3330 0 : for(size_t ip = 0; ip < numIP; ++ip)
3331 : {
3332 :
3333 : // compute approximated solution at integration point
3334 : number approxSolIP = 0.0;
3335 0 : for(size_t sh = 0; sh < num_sh; ++sh)
3336 : {
3337 : // get value at shape point (e.g. corner for P1 fct)
3338 : // and add shape fct at ip * value at shape
3339 0 : const number valSH = DoFRef((*m_pGridFct), ind[sh]);
3340 0 : approxSolIP += valSH * rTrialSpace.shape(sh, vLocIP[ip]);
3341 : }
3342 :
3343 : // get function value at ip
3344 0 : vValue[ip] = approxSolIP;
3345 :
3346 : }
3347 :
3348 : }
3349 0 : UG_CATCH_THROW("StdFuncIntegrand::evaluate: trial space missing.");
3350 0 : }
3351 : };
3352 :
3353 :
3354 : template <typename TGridFunction>
3355 0 : number StdFuncIntegralOnVertex(TGridFunction& gridFct,
3356 : size_t fct,
3357 : int si)
3358 : {
3359 : // integrate elements of subset
3360 : typedef typename TGridFunction::template dim_traits<0>::grid_base_object grid_base_object;
3361 : typedef typename TGridFunction::template dim_traits<0>::const_iterator const_iterator;
3362 :
3363 : // reset the result
3364 : number integral = 0;
3365 :
3366 : // note: this iterator is for the base elements, e.g. Face and not
3367 : // for the special type, e.g. Triangle, Quadrilateral
3368 : const_iterator iter = gridFct.template begin<grid_base_object>(si);
3369 : const_iterator iterEnd = gridFct.template end<grid_base_object>(si);
3370 :
3371 : // iterate over all elements
3372 0 : for(; iter != iterEnd; ++iter)
3373 : {
3374 : // get element
3375 : grid_base_object* pElem = *iter;
3376 :
3377 : std::vector<DoFIndex> ind; // aux. index array
3378 : gridFct.dof_indices(pElem, fct, ind);
3379 :
3380 : // compute approximated solution at integration point
3381 : number value = 0.0;
3382 0 : for(size_t sh = 0; sh < ind.size(); ++sh)
3383 : {
3384 0 : value += DoFRef(gridFct, ind[sh]);
3385 : }
3386 :
3387 : // add to global sum
3388 0 : integral += value;
3389 :
3390 : } // end elem
3391 :
3392 : // return the summed integral contributions of all elements
3393 0 : return integral;
3394 : }
3395 :
3396 : template <typename TGridFunction>
3397 : number StdFuncIntegralOnVertex(SmartPtr<TGridFunction> spGridFct, size_t fct, int si)
3398 : { return StdFuncIntegralOnVertex(*spGridFct, fct, si); }
3399 :
3400 :
3401 : template <typename TGridFunction>
3402 0 : number Integral(TGridFunction& gridFct, const char* cmp,
3403 : const char* subsets, int quadOrder)
3404 : {
3405 : // get function id of name
3406 : const size_t fct = gridFct.fct_id_by_name(cmp);
3407 :
3408 : // check that function exists
3409 0 : UG_COND_THROW(fct >= gridFct.num_fct(), "L2Norm: Function space does not contain"
3410 : " a function with name " << cmp << ".");
3411 :
3412 : // read subsets
3413 0 : SubsetGroup ssGrp(gridFct.domain()->subset_handler());
3414 0 : if(subsets != NULL)
3415 : {
3416 0 : ssGrp.add(TokenizeString(subsets));
3417 0 : if(!SameDimensionsInAllSubsets(ssGrp))
3418 0 : UG_THROW("IntegrateSubsets: Subsets '"<<subsets<<"' do not have same dimension."
3419 : "Can not integrate on subsets of different dimensions.");
3420 : }
3421 : else
3422 : {
3423 : // add all subsets and remove lower dim subsets afterwards
3424 0 : ssGrp.add_all();
3425 0 : RemoveLowerDimSubsets(ssGrp);
3426 : }
3427 :
3428 : // \TODO: This should be generalite in IntegrateSubset instead of doing it directly here
3429 : bool bOnlyVertex = true;
3430 0 : for(size_t s = 0; s < ssGrp.size(); ++s)
3431 0 : if(ssGrp.dim(s) != 0) bOnlyVertex = false;
3432 :
3433 0 : if(bOnlyVertex)
3434 : {
3435 : number value = 0;
3436 0 : for(size_t s = 0; s < ssGrp.size(); ++s)
3437 0 : value += StdFuncIntegralOnVertex(gridFct, fct, ssGrp[s]);
3438 :
3439 : #ifdef UG_PARALLEL
3440 : // sum over processes
3441 : if(pcl::NumProcs() > 1)
3442 : {
3443 : pcl::ProcessCommunicator com;
3444 : number local = value;
3445 : com.allreduce(&local, &value, 1, PCL_DT_DOUBLE, PCL_RO_SUM);
3446 : }
3447 : #endif
3448 0 : return value;
3449 : }
3450 :
3451 : StdFuncIntegrand<TGridFunction> integrand(&gridFct, fct);
3452 0 : return IntegrateSubsets(integrand, gridFct, subsets, quadOrder);
3453 : }
3454 :
3455 : template <typename TGridFunction>
3456 0 : number Integral(SmartPtr<TGridFunction> spGridFct, const char* cmp,
3457 : const char* subsets, int quadOrder)
3458 0 : { return Integral(*spGridFct, cmp, subsets, quadOrder); }
3459 :
3460 :
3461 : template <typename TGridFunction>
3462 0 : number Integral(SmartPtr<TGridFunction> spGridFct, const char* cmp,
3463 : const char* subsets)
3464 : {
3465 0 : return Integral(spGridFct, cmp, subsets, 1);
3466 : }
3467 : template <typename TGridFunction>
3468 0 : number Integral(SmartPtr<TGridFunction> spGridFct, const char* cmp)
3469 : {
3470 0 : return Integral(spGridFct, cmp, NULL, 1);
3471 : }
3472 :
3473 :
3474 : ////////////////////////////////////////////////////////////////////////////////
3475 : ////////////////////////////////////////////////////////////////////////////////
3476 : // Generic Boundary Integration Routine
3477 : ////////////////////////////////////////////////////////////////////////////////
3478 : ////////////////////////////////////////////////////////////////////////////////
3479 :
3480 : template <int WorldDim, int dim, typename TConstIterator>
3481 0 : number IntegralNormalComponentOnManifoldUsingFV1Geom(TConstIterator iterBegin,
3482 : TConstIterator iterEnd,
3483 : typename domain_traits<WorldDim>::position_accessor_type& aaPos,
3484 : const ISubsetHandler* ish,
3485 : IIntegrand<MathVector<WorldDim>, WorldDim>& integrand,
3486 : const SubsetGroup& bndSSGrp)
3487 : {
3488 : // reset the result
3489 : number integral = 0;
3490 :
3491 : // note: this iterator is for the base elements, e.g. Face and not
3492 : // for the special type, e.g. Triangle, Quadrilateral
3493 : TConstIterator iter = iterBegin;
3494 :
3495 : // this is the base element type (e.g. Face). This is the type when the
3496 : // iterators above are dereferenciated.
3497 : typedef typename domain_traits<dim>::grid_base_object grid_base_object;
3498 :
3499 : // create a FV1 Geometry
3500 0 : DimFV1Geometry<dim> geo;
3501 :
3502 : // specify, which subsets are boundary
3503 0 : for(size_t s = 0; s < bndSSGrp.size(); ++s)
3504 : {
3505 : // get subset index
3506 : const int bndSubset = bndSSGrp[s];
3507 :
3508 : // request this subset index as boundary subset. This will force the
3509 : // creation of boundary subsets when calling geo.update
3510 0 : geo.add_boundary_subset(bndSubset);
3511 : }
3512 :
3513 : // vector of corner coordinates of element corners (to be filled for each elem)
3514 : std::vector<MathVector<WorldDim> > vCorner;
3515 :
3516 : // iterate over all elements
3517 0 : for(; iter != iterEnd; ++iter)
3518 : {
3519 : // get element
3520 : grid_base_object* pElem = *iter;
3521 :
3522 : // get all corner coordinates
3523 : CollectCornerCoordinates(vCorner, *pElem, aaPos, true);
3524 :
3525 : // compute bf and grads at bip for element
3526 : try{
3527 0 : geo.update(pElem, &vCorner[0], ish);
3528 : }
3529 0 : UG_CATCH_THROW("IntegralNormalComponentOnManifold: "
3530 : "Cannot update Finite Volume Geometry.");
3531 :
3532 : // specify, which subsets are boundary
3533 0 : for(size_t s = 0; s < bndSSGrp.size(); ++s)
3534 : {
3535 : // get subset index
3536 : const int bndSubset = bndSSGrp[s];
3537 :
3538 : // get all bf of this subset
3539 : typedef typename DimFV1Geometry<dim>::BF BF;
3540 : const std::vector<BF>& vBF = geo.bf(bndSubset);
3541 :
3542 : // loop boundary faces
3543 0 : for(size_t b = 0; b < vBF.size(); ++b)
3544 : {
3545 : // get bf
3546 : const BF& bf = vBF[b];
3547 :
3548 : // value
3549 : MathVector<WorldDim> value;
3550 :
3551 : // jacobian
3552 : MathMatrix<dim, WorldDim> JT;
3553 0 : Inverse(JT, bf.JTInv());
3554 :
3555 : // compute integrand values at integration points
3556 : try
3557 : {
3558 0 : integrand.values(&value, &(bf.global_ip()),
3559 : pElem, &vCorner[0], &(bf.local_ip()),
3560 : &(JT),
3561 : 1);
3562 : }
3563 0 : UG_CATCH_THROW("IntegralNormalComponentOnManifold: Unable to compute values of "
3564 : "integrand at integration point.");
3565 :
3566 : // sum up contribution (normal includes area)
3567 0 : integral += VecDot(value, bf.normal());
3568 :
3569 : } // end bf
3570 : } // end bnd subsets
3571 : } // end elem
3572 :
3573 : // return the summed integral contributions of all elements
3574 0 : return integral;
3575 0 : }
3576 :
3577 :
3578 : template <int WorldDim, int dim, typename TConstIterator>
3579 0 : number IntegralNormalComponentOnManifoldGeneral(
3580 : TConstIterator iterBegin,
3581 : TConstIterator iterEnd,
3582 : typename domain_traits<WorldDim>::position_accessor_type& aaPos,
3583 : const ISubsetHandler* ish,
3584 : IIntegrand<MathVector<WorldDim>, WorldDim>& integrand,
3585 : const SubsetGroup& bndSSGrp,
3586 : int quadOrder,
3587 : Grid& grid)
3588 : {
3589 : // reset the result
3590 : number integral = 0;
3591 :
3592 : // note: this iterator is for the base elements, e.g. Face and not
3593 : // for the special type, e.g. Triangle, Quadrilateral
3594 : TConstIterator iter = iterBegin;
3595 :
3596 : // this is the base element type (e.g. Face). This is the type when the
3597 : // iterators above are dereferenciated.
3598 : typedef typename domain_traits<dim>::element_type Element;
3599 : typedef typename domain_traits<dim>::side_type Side;
3600 :
3601 : // vector of corner coordinates of element corners (to be filled for each elem)
3602 : std::vector<MathVector<WorldDim> > vCorner;
3603 : std::vector<int> vSubsetIndex;
3604 :
3605 : // iterate over all elements
3606 0 : for(; iter != iterEnd; ++iter)
3607 : {
3608 : // get element
3609 : Element* pElem = *iter;
3610 :
3611 : // get all corner coordinates
3612 : CollectCornerCoordinates(vCorner, *pElem, aaPos, true);
3613 :
3614 : // get reference object id
3615 0 : const ReferenceObjectID elemRoid = pElem->reference_object_id();
3616 :
3617 : // get sides
3618 : typename Grid::traits<Side>::secure_container vSide;
3619 : grid.associated_elements_sorted(vSide, pElem);
3620 0 : vSubsetIndex.resize(vSide.size());
3621 0 : for(size_t i = 0; i < vSide.size(); ++i)
3622 0 : vSubsetIndex[i] = ish->get_subset_index(vSide[i]);
3623 :
3624 : DimReferenceMapping<dim, WorldDim>& rMapping
3625 : = ReferenceMappingProvider::get<dim, WorldDim>(elemRoid, vCorner);
3626 :
3627 : const DimReferenceElement<dim>& rRefElem
3628 : = ReferenceElementProvider::get<dim>(elemRoid);
3629 :
3630 : // loop sub elements
3631 0 : for(size_t side = 0; side < vSide.size(); ++side)
3632 : {
3633 : // check if side used
3634 0 : if(!bndSSGrp.contains(vSubsetIndex[side])) continue;
3635 :
3636 : // get side
3637 : Side* pSide = vSide[side];
3638 :
3639 0 : std::vector<MathVector<WorldDim> > vSideCorner(rRefElem.num(dim-1, side, 0));
3640 0 : std::vector<MathVector<dim> > vLocalSideCorner(rRefElem.num(dim-1, side, 0));
3641 0 : for(size_t co = 0; co < vSideCorner.size(); ++co){
3642 0 : vSideCorner[co] = vCorner[rRefElem.id(dim-1, side, 0, co)];
3643 : vLocalSideCorner[co] = rRefElem.corner(rRefElem.id(dim-1, side, 0, co));
3644 : }
3645 :
3646 : // side quad rule
3647 0 : const ReferenceObjectID sideRoid = pSide->reference_object_id();
3648 : const QuadratureRule<dim-1>& rSideQuadRule
3649 0 : = QuadratureRuleProvider<dim-1>::get(sideRoid, quadOrder);
3650 :
3651 : // normal
3652 : MathVector<WorldDim> Normal;
3653 0 : ElementNormal<WorldDim>(sideRoid, Normal, &vSideCorner[0]);
3654 :
3655 : // quadrature points
3656 : const number* vWeight = rSideQuadRule.weights();
3657 : const size_t nip = rSideQuadRule.size();
3658 0 : std::vector<MathVector<dim> > vLocalIP(nip);
3659 0 : std::vector<MathVector<dim> > vGlobalIP(nip);
3660 :
3661 : DimReferenceMapping<dim-1, dim>& map
3662 : = ReferenceMappingProvider::get<dim-1, dim>(sideRoid, vLocalSideCorner);
3663 :
3664 0 : for(size_t ip = 0; ip < nip; ++ip)
3665 0 : map.local_to_global(vLocalIP[ip], rSideQuadRule.point(ip));
3666 :
3667 0 : for(size_t ip = 0; ip < nip; ++ip)
3668 0 : rMapping.local_to_global(vGlobalIP[ip], vLocalIP[ip]);
3669 :
3670 : // compute transformation matrices
3671 0 : std::vector<MathMatrix<dim-1, WorldDim> > vJT(nip);
3672 0 : map.jacobian_transposed(&(vJT[0]), rSideQuadRule.points(), nip);
3673 :
3674 0 : std::vector<MathMatrix<dim, WorldDim> > vElemJT(nip);
3675 0 : rMapping.jacobian_transposed(&(vElemJT[0]), &vLocalIP[0], nip);
3676 :
3677 0 : std::vector<MathVector<WorldDim> > vValue(nip);
3678 :
3679 : // compute integrand values at integration points
3680 : try
3681 : {
3682 0 : integrand.values(&vValue[0], &vGlobalIP[0],
3683 : pElem, &vCorner[0], &vLocalIP[0],
3684 : &(vElemJT[0]),
3685 : nip);
3686 : }
3687 0 : UG_CATCH_THROW("IntegralNormalComponentOnManifold: Unable to compute values of "
3688 : "integrand at integration point.");
3689 :
3690 : // loop integration points
3691 0 : for(size_t ip = 0; ip < nip; ++ip)
3692 : {
3693 : // get quadrature weight
3694 0 : const number weightIP = vWeight[ip];
3695 :
3696 : // get determinate of mapping
3697 0 : const number det = SqrtGramDeterminant(vJT[ip]);
3698 :
3699 : // add contribution of integration point
3700 0 : integral += VecDot(vValue[ip], Normal) * weightIP * det;
3701 : }
3702 : } // end bf
3703 : } // end elem
3704 :
3705 : // return the summed integral contributions of all elements
3706 0 : return integral;
3707 0 : }
3708 :
3709 : template <typename TGridFunction, int dim>
3710 0 : number IntegralNormalComponentOnManifoldSubset(SmartPtr<IIntegrand<MathVector<TGridFunction::dim>, TGridFunction::dim> > spIntegrand,
3711 : SmartPtr<TGridFunction> spGridFct,
3712 : int si, const SubsetGroup& bndSSGrp, int quadOrder)
3713 : {
3714 : // integrate elements of subset
3715 : typedef typename TGridFunction::template dim_traits<dim>::grid_base_object grid_base_object;
3716 : typedef typename TGridFunction::template dim_traits<dim>::const_iterator const_iterator;
3717 : static const int WorldDim = TGridFunction::dim;
3718 :
3719 0 : spIntegrand->set_subset(si);
3720 :
3721 0 : if(quadOrder == 1)
3722 : return IntegralNormalComponentOnManifoldUsingFV1Geom<WorldDim,dim,const_iterator>
3723 0 : (spGridFct->template begin<grid_base_object>(si),
3724 : spGridFct->template end<grid_base_object>(si),
3725 0 : spGridFct->domain()->position_accessor(),
3726 0 : spGridFct->domain()->subset_handler().get(),
3727 : *spIntegrand, bndSSGrp);
3728 : else{
3729 : UG_LOG(" #### IntegralNormalComponentOnManifoldSubset ####:\n")
3730 : return IntegralNormalComponentOnManifoldGeneral<WorldDim,dim,const_iterator>
3731 0 : (spGridFct->template begin<grid_base_object>(si),
3732 : spGridFct->template end<grid_base_object>(si),
3733 0 : spGridFct->domain()->position_accessor(),
3734 0 : spGridFct->domain()->subset_handler().get(),
3735 0 : *spIntegrand, bndSSGrp, quadOrder, *spGridFct->domain()->grid());
3736 : }
3737 : }
3738 :
3739 : template <typename TGridFunction>
3740 0 : number IntegralNormalComponentOnManifoldSubsets(
3741 : SmartPtr<IIntegrand<MathVector<TGridFunction::dim>, TGridFunction::dim> > spIntegrand,
3742 : SmartPtr<TGridFunction> spGridFct,
3743 : const char* BndSubsets, const char* InnerSubsets,
3744 : int quadOrder)
3745 : {
3746 : // world dimensions
3747 : static const int dim = TGridFunction::dim;
3748 :
3749 : // read subsets
3750 0 : SubsetGroup innerSSGrp(spGridFct->domain()->subset_handler());
3751 0 : if(InnerSubsets != NULL)
3752 : {
3753 0 : innerSSGrp.add(TokenizeString(InnerSubsets));
3754 0 : if(!SameDimensionsInAllSubsets(innerSSGrp))
3755 0 : UG_THROW("IntegralNormalComponentOnManifold: Subsets '"<<InnerSubsets<<"' do not have same dimension."
3756 : "Can not integrate on subsets of different dimensions.");
3757 : }
3758 : else
3759 : {
3760 : // add all subsets and remove lower dim subsets afterwards
3761 0 : innerSSGrp.add_all();
3762 0 : RemoveLowerDimSubsets(innerSSGrp);
3763 : }
3764 :
3765 : // read subsets
3766 0 : SubsetGroup bndSSGrp(spGridFct->domain()->subset_handler());
3767 0 : if(BndSubsets != NULL)
3768 0 : bndSSGrp.add(TokenizeString(BndSubsets));
3769 : else
3770 0 : UG_THROW("IntegralNormalComponentOnManifold: No boundary subsets passed.");
3771 :
3772 : // reset value
3773 : number value = 0;
3774 :
3775 : // loop subsets
3776 0 : for(size_t i = 0; i < innerSSGrp.size(); ++i)
3777 : {
3778 : // get subset index
3779 : const int si = innerSSGrp[i];
3780 :
3781 : // check dimension
3782 0 : if(innerSSGrp.dim(i) != dim && innerSSGrp.dim(i) != DIM_SUBSET_EMPTY_GRID)
3783 0 : UG_THROW("IntegralNormalComponentOnManifold: Dimension of inner subset is "<<
3784 : innerSSGrp.dim(i)<<", but only World Dimension "<<dim<<
3785 : " subsets can be used for inner subsets.");
3786 :
3787 : // integrate elements of subset
3788 0 : switch(innerSSGrp.dim(i))
3789 : {
3790 : case DIM_SUBSET_EMPTY_GRID: break;
3791 0 : case dim: value += IntegralNormalComponentOnManifoldSubset<TGridFunction, dim>(spIntegrand, spGridFct, si, bndSSGrp, quadOrder); break;
3792 0 : default: UG_THROW("IntegralNormalComponentOnManifold: Dimension "<<innerSSGrp.dim(i)<<" not supported. "
3793 : " World dimension is "<<dim<<".");
3794 : }
3795 : }
3796 :
3797 : #ifdef UG_PARALLEL
3798 : // sum over processes
3799 : if(pcl::NumProcs() > 1)
3800 : {
3801 : pcl::ProcessCommunicator com;
3802 : number local = value;
3803 : com.allreduce(&local, &value, 1, PCL_DT_DOUBLE, PCL_RO_SUM);
3804 : }
3805 : #endif
3806 :
3807 : // return the result
3808 0 : return value;
3809 : }
3810 :
3811 : template <typename TGridFunction>
3812 0 : number IntegralNormalComponentOnManifold(
3813 : SmartPtr<UserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > spData,
3814 : SmartPtr<TGridFunction> spGridFct,
3815 : const char* BndSubset, const char* InnerSubset,
3816 : number time,
3817 : int quadOrder)
3818 : {
3819 0 : SmartPtr<IIntegrand<MathVector<TGridFunction::dim>, TGridFunction::dim> > spIntegrand
3820 0 : = make_sp(new UserDataIntegrand<MathVector<TGridFunction::dim>, TGridFunction>(spData, &(*spGridFct), time));
3821 :
3822 0 : return IntegralNormalComponentOnManifoldSubsets(spIntegrand, spGridFct, BndSubset, InnerSubset, quadOrder);
3823 : }
3824 :
3825 : template <typename TGridFunction>
3826 0 : number IntegralNormalComponentOnManifold(
3827 : SmartPtr<UserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > spData,
3828 : SmartPtr<TGridFunction> spGridFct,
3829 : const char* BndSubset, const char* InnerSubset,
3830 : number time)
3831 0 : {return IntegralNormalComponentOnManifold(spData, spGridFct, BndSubset, InnerSubset, time, 1);}
3832 :
3833 : template <typename TGridFunction>
3834 0 : number IntegralNormalComponentOnManifold(
3835 : SmartPtr<UserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > spData,
3836 : SmartPtr<TGridFunction> spGridFct,
3837 : const char* BndSubset,
3838 : number time)
3839 0 : {return IntegralNormalComponentOnManifold(spData, spGridFct, BndSubset, NULL, time, 1);}
3840 :
3841 : template <typename TGridFunction>
3842 0 : number IntegralNormalComponentOnManifold(
3843 : SmartPtr<UserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > spData,
3844 : SmartPtr<TGridFunction> spGridFct,
3845 : const char* BndSubset, const char* InnerSubset)
3846 0 : {return IntegralNormalComponentOnManifold(spData, spGridFct, BndSubset, InnerSubset, 0.0, 1);}
3847 :
3848 : template <typename TGridFunction>
3849 0 : number IntegralNormalComponentOnManifold(
3850 : SmartPtr<UserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > spData,
3851 : SmartPtr<TGridFunction> spGridFct,
3852 : const char* BndSubset)
3853 0 : {return IntegralNormalComponentOnManifold(spData, spGridFct, BndSubset, NULL, 0.0, 1);}
3854 :
3855 : ///////////////
3856 : // lua data
3857 : ///////////////
3858 :
3859 : #ifdef UG_FOR_LUA
3860 : template <typename TGridFunction>
3861 0 : number IntegralNormalComponentOnManifold(
3862 : const char* luaFct, SmartPtr<TGridFunction> spGridFct,
3863 : const char* BndSubset, const char* InnerSubset,
3864 : number time, int quadOrder)
3865 : {
3866 : static const int dim = TGridFunction::dim;
3867 0 : SmartPtr<UserData<MathVector<dim>, dim> > sp =
3868 : LuaUserDataFactory<MathVector<dim>, dim>::create(luaFct);
3869 0 : return IntegralNormalComponentOnManifold(sp, spGridFct, BndSubset, InnerSubset, time, quadOrder);
3870 : }
3871 :
3872 : template <typename TGridFunction>
3873 0 : number IntegralNormalComponentOnManifold(
3874 : const char* luaFct, SmartPtr<TGridFunction> spGridFct,
3875 : const char* BndSubset, const char* InnerSubset,
3876 : number time)
3877 0 : {return IntegralNormalComponentOnManifold(luaFct, spGridFct, BndSubset, InnerSubset, time, 1);}
3878 :
3879 : template <typename TGridFunction>
3880 0 : number IntegralNormalComponentOnManifold(
3881 : const char* luaFct, SmartPtr<TGridFunction> spGridFct,
3882 : const char* BndSubset,
3883 : number time)
3884 0 : {return IntegralNormalComponentOnManifold(luaFct, spGridFct, BndSubset, NULL, time, 1);}
3885 :
3886 : template <typename TGridFunction>
3887 0 : number IntegralNormalComponentOnManifold(
3888 : const char* luaFct, SmartPtr<TGridFunction> spGridFct,
3889 : const char* BndSubset, const char* InnerSubset)
3890 0 : {return IntegralNormalComponentOnManifold(luaFct, spGridFct, BndSubset, InnerSubset, 0.0, 1);}
3891 :
3892 : template <typename TGridFunction>
3893 0 : number IntegralNormalComponentOnManifold(
3894 : const char* luaFct, SmartPtr<TGridFunction> spGridFct,
3895 : const char* BndSubset)
3896 0 : {return IntegralNormalComponentOnManifold(luaFct, spGridFct, BndSubset, NULL, 0.0, 1);}
3897 : #endif
3898 :
3899 : ////////////////////////////////////////////////////////////////////////////////
3900 : // Boundary Integral
3901 : ////////////////////////////////////////////////////////////////////////////////
3902 : /// Integrates a component over some boundary subsets
3903 : /**
3904 : * This function integrates \f$ -c \cdot \vec{n} \f$ over some selected
3905 : * boundary subsets, where c is a component of the grid function u. The
3906 : * integral is taken over all boundary subset manifold geometric objects
3907 : *
3908 : * @param u a grid function
3909 : * @param cmp the component c which should be integrated
3910 : * @param BndSubset a comma-separated string of symbolic boundary subset names
3911 : * @return the integral
3912 : */
3913 :
3914 : template <typename TGridFunction>
3915 0 : number IntegrateNormalComponentOnManifold(TGridFunction& u, const char* cmp,
3916 : const char* BndSubset)
3917 : {
3918 : // get function id of name
3919 : const size_t fct = u.fct_id_by_name(cmp);
3920 :
3921 : // check that function exists
3922 0 : if(fct >= u.num_fct())
3923 0 : UG_THROW("IntegrateNormalComponentOnManifold: Function space does not contain"
3924 : " a function with name " << cmp << ".");
3925 :
3926 : if(u.local_finite_element_id(fct) != LFEID(LFEID::LAGRANGE, TGridFunction::dim, 1))
3927 0 : UG_THROW("IntegrateNormalComponentOnManifold:"
3928 : "Only implemented for Lagrange P1 functions.");
3929 :
3930 : // read bnd subsets
3931 0 : SubsetGroup bndSSGrp(u.domain()->subset_handler());
3932 0 : bndSSGrp.add(TokenizeString(BndSubset));
3933 :
3934 : // reset value
3935 : number value = 0;
3936 :
3937 : // loop subsets
3938 0 : for(size_t i = 0; i < bndSSGrp.size(); ++i)
3939 : {
3940 : // get subset index
3941 : const int si = bndSSGrp[i];
3942 :
3943 : // skip if function is not defined in subset
3944 0 : if(!u.is_def_in_subset(fct, si)) continue;
3945 :
3946 : // create integration kernel
3947 : static const int dim = TGridFunction::dim;
3948 :
3949 : // integrate elements of subset
3950 : typedef typename domain_traits<dim>::grid_base_object grid_base_object;
3951 :
3952 : // get iterators for all elems on subset
3953 : typedef typename TGridFunction::template dim_traits<dim>::const_iterator const_iterator;
3954 : const_iterator iter = u.template begin<grid_base_object>();
3955 : const_iterator iterEnd = u.template end<grid_base_object>();
3956 :
3957 : // create a FV1 Geometry
3958 0 : DimFV1Geometry<dim> geo;
3959 :
3960 : // specify, which subsets are boundary
3961 0 : for(size_t s = 0; s < bndSSGrp.size(); ++s)
3962 : {
3963 : // get subset index
3964 : const int bndSubset = bndSSGrp[s];
3965 :
3966 : // request this subset index as boundary subset. This will force the
3967 : // creation of boundary subsets when calling geo.update
3968 0 : geo.add_boundary_subset(bndSubset);
3969 : }
3970 :
3971 : // vector of corner coordinates of element corners (to be filled for each elem)
3972 : std::vector<MathVector<dim> > vCorner;
3973 :
3974 : // loop elements of subset
3975 0 : for( ; iter != iterEnd; ++iter)
3976 : {
3977 : // get element
3978 : grid_base_object* elem = *iter;
3979 :
3980 : // get all corner coordinates
3981 0 : CollectCornerCoordinates(vCorner, *elem, u.domain()->position_accessor(), true);
3982 :
3983 : // compute bf and grads at bip for element
3984 : try{
3985 0 : geo.update(elem, &vCorner[0], u.domain()->subset_handler().get());
3986 : }
3987 0 : UG_CATCH_THROW("IntegrateNormalComponenOnManifold: "
3988 : "Cannot update Finite Volume Geometry.");
3989 :
3990 : // get fct multi-indices of element
3991 : std::vector<DoFIndex> ind;
3992 : u.dof_indices(elem, fct, ind);
3993 :
3994 : // specify, which subsets are boundary
3995 0 : for(size_t s = 0; s < bndSSGrp.size(); ++s)
3996 : {
3997 : // get subset index
3998 : const int bndSubset = bndSSGrp[s];
3999 :
4000 : // get all bf of this subset
4001 : typedef typename DimFV1Geometry<dim>::BF BF;
4002 : const std::vector<BF>& vBF = geo.bf(bndSubset);
4003 :
4004 : // loop boundary faces
4005 0 : for(size_t b = 0; b < vBF.size(); ++b)
4006 : {
4007 : // get bf
4008 : const BF& bf = vBF[b];
4009 :
4010 : // get normal on bf
4011 : const MathVector<dim>& normal = bf.normal();
4012 :
4013 : // check multi indices
4014 : UG_ASSERT(ind.size() == bf.num_sh(),
4015 : "IntegrateNormalComponentOnManifold::values: Wrong number of"
4016 : " multi indices, ind: "<<ind.size() << ", bf.num_sh: "
4017 : << bf.num_sh());
4018 :
4019 : MathVector<dim> val; VecSet(val, 0.0);
4020 0 : for(size_t sh = 0; sh < bf.num_sh(); ++sh)
4021 : {
4022 0 : const number fctVal = DoFRef(u, ind[sh]);
4023 : const MathVector<dim>& ip = bf.global_ip();
4024 : VecScaleAdd(val, 1.0, val, fctVal, ip);
4025 : }
4026 :
4027 : // sum up contributions
4028 0 : value -= VecDot(val, normal);
4029 : }
4030 : }
4031 : }
4032 : }
4033 :
4034 : #ifdef UG_PARALLEL
4035 : // sum over processes
4036 : if(pcl::NumProcs() > 1)
4037 : {
4038 : pcl::ProcessCommunicator com;
4039 : number local = value;
4040 : com.allreduce(&local, &value, 1, PCL_DT_DOUBLE, PCL_RO_SUM);
4041 : }
4042 : #endif
4043 :
4044 : // return the result
4045 0 : return value;
4046 : }
4047 :
4048 : /// Integrates the Flux of a component over some boundary subsets
4049 : /**
4050 : * This function integrates \f$ - \nabla c \cdot \vec{n} \f$ over some selected
4051 : * boundary subsets. In order to compute the gradient of a function a world-
4052 : * dimension element must be given and, thus, some "inner" subsets have to be
4053 : * specified to indicate the full-dimensional subsets. The integral sum is then
4054 : * taken over all boundary subset manifold geometric objects, that are part of
4055 : * an element of the scheduled "inner" elements
4056 : *
4057 : * @param u a grid function
4058 : * @param cmp the component, whose gradient should be integrated
4059 : * @param BndSubset a comma-separated string of symbolic boundary subset names
4060 : * @param InnerSubset a comma-separated string of symbolic inner subset names
4061 : * @return the integral
4062 : */
4063 : template <typename TGridFunction>
4064 0 : number IntegrateNormalGradientOnManifold(TGridFunction& u, const char* cmp,
4065 : const char* BndSubset, const char* InnerSubset = NULL)
4066 : {
4067 : // get function id of name
4068 : const size_t fct = u.fct_id_by_name(cmp);
4069 :
4070 : // check that function exists
4071 0 : if(fct >= u.num_fct())
4072 0 : UG_THROW("IntegrateNormalGradientOnManifold: Function space does not contain"
4073 : " a function with name " << cmp << ".");
4074 :
4075 : if(u.local_finite_element_id(fct) != LFEID(LFEID::LAGRANGE, TGridFunction::dim, 1))
4076 0 : UG_THROW("IntegrateNormalGradientOnManifold:"
4077 : "Only implemented for Lagrange P1 functions.");
4078 :
4079 : // read subsets
4080 0 : SubsetGroup innerSSGrp(u.domain()->subset_handler());
4081 0 : if(InnerSubset != NULL)
4082 0 : innerSSGrp.add(TokenizeString(InnerSubset));
4083 : else // add all if no subset specified
4084 0 : innerSSGrp.add_all();
4085 :
4086 : // read bnd subsets
4087 0 : SubsetGroup bndSSGrp(u.domain()->subset_handler());
4088 0 : if(InnerSubset != NULL){
4089 0 : bndSSGrp.add(TokenizeString(BndSubset));
4090 : }
4091 : else{
4092 0 : UG_THROW("IntegrateNormalGradientOnManifold: No boundary subsets specified. Aborting.");
4093 : }
4094 :
4095 : // reset value
4096 : number value = 0;
4097 :
4098 : // loop subsets
4099 0 : for(size_t i = 0; i < innerSSGrp.size(); ++i)
4100 : {
4101 : // get subset index
4102 : const int si = innerSSGrp[i];
4103 :
4104 : // skip if function is not defined in subset
4105 0 : if(!u.is_def_in_subset(fct, si)) continue;
4106 :
4107 :
4108 0 : if (innerSSGrp.dim(i) != TGridFunction::dim)
4109 0 : UG_THROW("IntegrateNormalGradientOnManifold: Element dimension does not match world dimension!");
4110 :
4111 : // create integration kernel
4112 : static const int dim = TGridFunction::dim;
4113 :
4114 : // integrate elements of subset
4115 : typedef typename domain_traits<dim>::grid_base_object grid_base_object;
4116 :
4117 : // get iterators for all elems on subset
4118 : typedef typename TGridFunction::template dim_traits<dim>::const_iterator const_iterator;
4119 : const_iterator iter = u.template begin<grid_base_object>();
4120 : const_iterator iterEnd = u.template end<grid_base_object>();
4121 :
4122 : // create a FV1 Geometry
4123 0 : DimFV1Geometry<dim> geo;
4124 :
4125 : // specify, which subsets are boundary
4126 0 : for(size_t s = 0; s < bndSSGrp.size(); ++s)
4127 : {
4128 : // get subset index
4129 : const int bndSubset = bndSSGrp[s];
4130 :
4131 : // request this subset index as boundary subset. This will force the
4132 : // creation of boundary subsets when calling geo.update
4133 0 : geo.add_boundary_subset(bndSubset);
4134 : }
4135 :
4136 : // vector of corner coordinates of element corners (to be filled for each elem)
4137 : std::vector<MathVector<dim> > vCorner;
4138 :
4139 : // loop elements of subset
4140 0 : for( ; iter != iterEnd; ++iter)
4141 : {
4142 : // get element
4143 : grid_base_object* elem = *iter;
4144 :
4145 : // get all corner coordinates
4146 0 : CollectCornerCoordinates(vCorner, *elem, u.domain()->position_accessor(), true);
4147 :
4148 : // compute bf and grads at bip for element
4149 : try{
4150 0 : geo.update(elem, &vCorner[0], u.domain()->subset_handler().get());
4151 : }
4152 0 : UG_CATCH_THROW("IntegrateNormalGradientOnManifold: "
4153 : "Cannot update Finite Volume Geometry.");
4154 :
4155 : // get fct multi-indeces of element
4156 : std::vector<DoFIndex> ind;
4157 : u.dof_indices(elem, fct, ind);
4158 :
4159 : // specify, which subsets are boundary
4160 0 : for(size_t s = 0; s < bndSSGrp.size(); ++s)
4161 : {
4162 : // get subset index
4163 : const int bndSubset = bndSSGrp[s];
4164 :
4165 : // get all bf of this subset
4166 : typedef typename DimFV1Geometry<dim>::BF BF;
4167 : const std::vector<BF>& vBF = geo.bf(bndSubset);
4168 :
4169 : // loop boundary faces
4170 0 : for(size_t b = 0; b < vBF.size(); ++b)
4171 : {
4172 : // get bf
4173 : const BF& bf = vBF[b];
4174 :
4175 : // get normal on bf
4176 : const MathVector<dim>& normal = bf.normal();
4177 :
4178 : // check multi indices
4179 : UG_ASSERT(ind.size() == bf.num_sh(),
4180 : "IntegrateNormalGradientOnManifold::values: Wrong number of"
4181 : " multi indices, ind: "<<ind.size() << ", bf.num_sh: "
4182 : << bf.num_sh());
4183 :
4184 : // compute gradient of solution at bip
4185 : MathVector<dim> grad; VecSet(grad, 0.0);
4186 0 : for(size_t sh = 0; sh < bf.num_sh(); ++sh)
4187 : {
4188 0 : const number fctVal = DoFRef(u, ind[sh]);
4189 :
4190 : VecScaleAdd(grad, 1.0, grad, fctVal, bf.global_grad(sh));
4191 : }
4192 :
4193 : // sum up contributions
4194 0 : value -= VecDot(grad, normal);
4195 : }
4196 : }
4197 : }
4198 : }
4199 :
4200 : #ifdef UG_PARALLEL
4201 : // sum over processes
4202 : if(pcl::NumProcs() > 1)
4203 : {
4204 : pcl::ProcessCommunicator com;
4205 : number local = value;
4206 : com.allreduce(&local, &value, 1, PCL_DT_DOUBLE, PCL_RO_SUM);
4207 : }
4208 : #endif
4209 :
4210 : // return the result
4211 0 : return value;
4212 : }
4213 :
4214 : } // namespace ug
4215 :
4216 : #endif /*__H__UG__LIB_DISC__FUNCTION_SPACES__INTEGRATE__*/
|
