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