Line data Source code
1 : /*
2 : * Copyright (c) 2012-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 §7):
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_SPACE__GRID_FUNCTION_USER_DATA__
34 : #define __H__UG__LIB_DISC__FUNCTION_SPACE__GRID_FUNCTION_USER_DATA__
35 :
36 : #include "common/common.h"
37 :
38 : #include "lib_grid/tools/subset_group.h"
39 :
40 : #include "lib_disc/common/function_group.h"
41 : #include "lib_disc/common/groups_util.h"
42 : #include "lib_disc/quadrature/quadrature.h"
43 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
44 : #include "lib_disc/spatial_disc/user_data/std_user_data.h"
45 : #include "lib_disc/reference_element/reference_mapping_provider.h"
46 :
47 :
48 : namespace ug{
49 :
50 : /**
51 : * Class for StdDependentUserData interface to scalar values from GridFunction objects.
52 : *
53 : * This class provides a StdDependentUserData interface to scalar values from GridFunction
54 : * objects so that they can be used in import parameters. Note that the GridFunction object
55 : * must have the same FunctionPattern as the owner of the import parameter.
56 : *
57 : * REMARK:
58 : * This class is used to reuse the unknowns of variables of discretizations. It provides
59 : * "derivatives w.r.t. the unknowns", too. Note that in this case, the GridFunction object
60 : * must be the same as the solution. An attempt to use a grid function based on a different
61 : * ApproximationSpace would lead to errors. Consider ExplicitGridFunctionValue
62 : * for this purpose instead.
63 : */
64 : template <typename TGridFunction>
65 : class GridFunctionNumberData
66 : : public StdDependentUserData<GridFunctionNumberData<TGridFunction>,
67 : number, TGridFunction::dim>
68 : {
69 : public:
70 : // world dimension of grid function
71 : static const int dim = TGridFunction::dim;
72 :
73 : private:
74 : // grid function
75 : SmartPtr<TGridFunction> m_spGridFct;
76 :
77 : // component of function
78 : size_t m_fct;
79 :
80 : // local finite element id
81 : LFEID m_lfeID;
82 :
83 : public:
84 : /// constructor
85 0 : GridFunctionNumberData(SmartPtr<TGridFunction> spGridFct, const char* cmp)
86 0 : {
87 0 : assign(spGridFct, cmp);
88 0 : }
89 :
90 : GridFunctionNumberData(){}
91 :
92 0 : void assign(SmartPtr<TGridFunction> spGridFct, const char* cmp)
93 : {
94 0 : m_spGridFct = spGridFct;
95 0 : this->set_functions(cmp);
96 :
97 : // get function id of name
98 0 : m_fct = spGridFct->fct_id_by_name(cmp);
99 :
100 : // check that function exists
101 0 : if(m_fct >= spGridFct->num_fct())
102 0 : UG_THROW("GridFunctionNumberData: Function space does not contain"
103 : " a function with name " << cmp << ".");
104 :
105 : // local finite element id
106 0 : m_lfeID = spGridFct->local_finite_element_id(m_fct);
107 0 : }
108 :
109 0 : virtual bool continuous() const
110 : {
111 0 : return LocalFiniteElementProvider::continuous(m_lfeID);
112 : }
113 :
114 : template <int refDim>
115 0 : void eval_and_deriv(number vValue[],
116 : const MathVector<dim> vGlobIP[],
117 : number time, int si,
118 : GridObject* elem,
119 : const MathVector<dim> vCornerCoords[],
120 : const MathVector<refDim> vLocIP[],
121 : const size_t nip,
122 : LocalVector* u,
123 : bool bDeriv,
124 : int s,
125 : std::vector<std::vector<number> > vvvDeriv[],
126 : const MathMatrix<refDim, dim>* vJT = NULL) const
127 : {
128 : // reference object id
129 0 : const ReferenceObjectID roid = elem->reference_object_id();
130 :
131 : // get trial space
132 : try{
133 : const LocalShapeFunctionSet<refDim>& rTrialSpace =
134 0 : LocalFiniteElementProvider::get<refDim>(roid, m_lfeID);
135 :
136 : // memory for shapes
137 : std::vector<number> vShape;
138 : // memory for indices
139 : std::vector<DoFIndex> ind;
140 :
141 : // loop ips
142 0 : for(size_t ip = 0; ip < nip; ++ip)
143 : {
144 : // evaluate at shapes at ip
145 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
146 :
147 : // get multiindices of element
148 0 : m_spGridFct->dof_indices(elem, m_fct, ind);
149 :
150 : // compute solution at integration point
151 0 : vValue[ip] = 0.0;
152 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
153 : {
154 0 : const number valSH = DoFRef(*m_spGridFct, ind[sh]);
155 0 : vValue[ip] += valSH * vShape[sh];
156 : }
157 : }
158 :
159 0 : if(bDeriv){
160 0 : for(size_t ip = 0; ip < nip; ++ip){
161 : // evaluate at shapes at ip
162 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
163 :
164 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
165 0 : vvvDeriv[ip][0][sh] = vShape[sh];
166 : }
167 : }
168 0 : }
169 0 : UG_CATCH_THROW("GridFunctionNumberData: Shape Function Set missing for"
170 : " Reference Object: "<<roid<<", Trial Space: "
171 : <<m_lfeID<<", refDim="<<refDim);
172 :
173 0 : }
174 :
175 0 : virtual void operator() (number& value,
176 : const MathVector<dim>& globIP,
177 : number time, int si,
178 : Vertex* vrt) const
179 : {
180 : std::vector<DoFIndex> ind;
181 0 : if (m_spGridFct->dof_indices(vrt, m_fct, ind) != 1)
182 0 : UG_THROW ("GridFunctionNumberData: Values at vertices of the grid function are not uniquely defined.");
183 0 : value = DoFRef(*m_spGridFct, ind[0]);
184 0 : }
185 : };
186 :
187 : /**
188 : * Class for StdDependentUserData interface to vectors from GridFunction objects.
189 : *
190 : * This class provides a StdDependentUserData interface to vectors from GridFunction
191 : * objects so that they can be used in import parameters. Note that the GridFunction object
192 : * must have the same FunctionPattern as the owner of the import parameter.
193 : *
194 : * REMARK:
195 : * This class is used to reuse the unknowns of variables of discretizations. It provides
196 : * "derivatives w.r.t. the unknowns", too. Note that in this case the GridFunction object
197 : * must be the same as the solution. An attempt to use a grid function based on a different
198 : * ApproximationSpace would lead to errors. Consider ExplicitGridFunctionVector
199 : * for this purpose instead.
200 : */
201 : template <typename TGridFunction>
202 : class GridFunctionVectorData
203 : : public StdDependentUserData<GridFunctionVectorData<TGridFunction>,
204 : MathVector<TGridFunction::dim>, TGridFunction::dim>
205 : {
206 : public:
207 : // world dimension of grid function
208 : static const int dim = TGridFunction::dim;
209 :
210 : private:
211 : // grid function
212 : SmartPtr<TGridFunction> m_spGridFct;
213 :
214 : // component of function
215 : size_t m_vfct[dim];
216 :
217 : // local finite element id
218 : LFEID m_vlfeID[dim];
219 :
220 : public:
221 : /// constructor
222 0 : GridFunctionVectorData(SmartPtr<TGridFunction> spGridFct, const char* cmp)
223 0 : : m_spGridFct(spGridFct)
224 : {
225 0 : this->set_functions(cmp);
226 :
227 : // create function group of this elem disc
228 : try{
229 : // get strings
230 0 : std::string fctString = std::string(cmp);
231 :
232 : // tokenize strings and select functions
233 : std::vector<std::string> tokens;
234 0 : TokenizeString(fctString, tokens, ',');
235 :
236 0 : for(size_t i = 0; i < tokens.size(); ++i)
237 0 : RemoveWhitespaceFromString(tokens[i]);
238 :
239 0 : if((int)tokens.size() != dim)
240 0 : UG_THROW("GridFunctionVectorData: Needed "<<dim<<" components "
241 : "in symbolic function names, but given: "<<cmp);
242 :
243 : // get function id of name
244 0 : for(int i = 0; i < dim; ++i){
245 0 : m_vfct[i] = spGridFct->fct_id_by_name(tokens[i].c_str());
246 0 : m_vlfeID[i] = spGridFct->local_finite_element_id(m_vfct[i]);
247 : }
248 :
249 0 : }UG_CATCH_THROW("GridFunctionVectorData: Cannot find some symbolic function "
250 : "name in '"<<cmp<<"'.");
251 0 : };
252 :
253 : GridFunctionVectorData(SmartPtr<TGridFunction> spGridFct, std::vector<std::string> vCmp)
254 : : m_spGridFct(spGridFct)
255 : {
256 : this->set_functions(vCmp);
257 :
258 : // create function group of this elem disc
259 : try
260 : {
261 : if ((int)vCmp.size() != dim)
262 : UG_THROW("GridFunctionVectorData: Needed "<<dim<<" components "
263 : "in symbolic function names, but given: "<< (int) vCmp.size());
264 :
265 : // get function id of name
266 : for (size_t i = 0; i < (size_t) dim; ++i)
267 : {
268 : m_vfct[i] = spGridFct->fct_id_by_name(vCmp[i].c_str());
269 : m_vlfeID[i] = spGridFct->local_finite_element_id(m_vfct[i]);
270 : }
271 :
272 : }UG_CATCH_THROW("GridFunctionVectorData: Cannot find some symbolic function "
273 : "name in component vector.");
274 : }
275 :
276 0 : virtual bool continuous() const
277 : {
278 0 : for(int i = 0; i < dim; ++i)
279 0 : if(!LocalFiniteElementProvider::continuous(m_vlfeID[i]))
280 : return false;
281 : return true;
282 : }
283 :
284 : template <int refDim>
285 0 : void eval_and_deriv(MathVector<dim> vValue[],
286 : const MathVector<dim> vGlobIP[],
287 : number time, int si,
288 : GridObject* elem,
289 : const MathVector<dim> vCornerCoords[],
290 : const MathVector<refDim> vLocIP[],
291 : const size_t nip,
292 : LocalVector* u,
293 : bool bDeriv,
294 : int s,
295 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[],
296 : const MathMatrix<refDim, dim>* vJT = NULL) const
297 : {
298 : // reference object id
299 0 : const ReferenceObjectID roid = elem->reference_object_id();
300 :
301 : // memory for shapes
302 : std::vector<number> vShape;
303 : // memory for indices
304 : std::vector<DoFIndex> ind;
305 :
306 : // loop ips
307 : try{
308 0 : for(size_t ip = 0; ip < nip; ++ip)
309 : {
310 0 : for(int d = 0; d < dim; ++d)
311 : {
312 : const LocalShapeFunctionSet<refDim>& rTrialSpace =
313 0 : LocalFiniteElementProvider::get<refDim>(roid, m_vlfeID[d]);
314 :
315 : // evaluate at shapes at ip
316 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
317 :
318 : // get multiindices of element
319 0 : m_spGridFct->dof_indices(elem, m_vfct[d], ind);
320 :
321 : // compute solution at integration point
322 0 : vValue[ip][d] = 0.0;
323 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
324 : {
325 0 : const number valSH = DoFRef( *m_spGridFct, ind[sh]);
326 0 : vValue[ip][d] += valSH * vShape[sh];
327 : }
328 : }
329 : }
330 :
331 :
332 0 : if(bDeriv){
333 0 : for(size_t ip = 0; ip < nip; ++ip){
334 0 : for(int d = 0; d < dim; ++d)
335 : {
336 : const LocalShapeFunctionSet<refDim>& rTrialSpace =
337 0 : LocalFiniteElementProvider::get<refDim>(roid, m_vlfeID[d]);
338 : // evaluate at shapes at ip
339 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
340 :
341 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
342 0 : vvvDeriv[ip][d][sh][d] = vShape[sh];
343 : }
344 : }
345 : }
346 :
347 :
348 0 : }UG_CATCH_THROW("GridFunctionVectorData: Reference Object: "
349 : <<roid<<", refDim="<<refDim);
350 0 : }
351 :
352 0 : virtual void operator() (MathVector<dim>& value,
353 : const MathVector<dim>& globIP,
354 : number time, int si,
355 : Vertex* vrt) const
356 : {
357 : std::vector<DoFIndex> ind;
358 0 : for(int d = 0; d < dim; ++d)
359 : {
360 0 : if (m_spGridFct->dof_indices(vrt, m_vfct[d], ind) != 1)
361 0 : UG_THROW ("GridFunctionVectorData: Values at vertices of the grid function are not uniquely defined.");
362 0 : value[d] = DoFRef(*m_spGridFct, ind[0]);
363 : }
364 0 : }
365 : };
366 :
367 :
368 : /**
369 : * Class for StdDependentUserData interface to the gradient of scalar values from GridFunction objects.
370 : *
371 : * This class provides a StdDependentUserData interface to the gradient of scalar values
372 : * from GridFunction objects so that they can be used in import parameters. Note that the
373 : * GridFunction object must have the same FunctionPattern as the owner of the import parameter.
374 : *
375 : * REMARK:
376 : * This class is used to reuse the unknowns of variables of discretizations. It provides
377 : * "derivatives w.r.t. the unknowns", too. Note that in this case the GridFunction object
378 : * must be the same as the solution. An attempt to use a grid function based on a different
379 : * ApproximationSpace would lead to errors. Consider ExplicitGridFunctionGradient
380 : * for this purpose instead.
381 : */
382 : template <typename TGridFunction>
383 : class GridFunctionGradientData
384 : : public StdDependentUserData<GridFunctionGradientData<TGridFunction> ,
385 : MathVector<TGridFunction::dim>, TGridFunction::dim>
386 : {
387 : public:
388 : // world dimension of grid function
389 : static const int dim = TGridFunction::dim;
390 :
391 : private:
392 : // grid function
393 : SmartPtr<TGridFunction> m_spGridFct;
394 :
395 : // component of function
396 : size_t m_fct;
397 :
398 : // local finite element id
399 : LFEID m_lfeID;
400 :
401 : public:
402 : /// constructor
403 0 : GridFunctionGradientData(SmartPtr<TGridFunction> spGridFct, const char* cmp)
404 0 : : m_spGridFct(spGridFct)
405 : {
406 0 : this->set_functions(cmp);
407 :
408 : // get function id of name
409 0 : m_fct = spGridFct->fct_id_by_name(cmp);
410 :
411 : // check that function exists
412 0 : if(m_fct >= spGridFct->num_fct())
413 0 : UG_THROW("GridFunctionGradientData: Function space does not contain"
414 : " a function with name " << cmp << ".");
415 :
416 : // local finite element id
417 0 : m_lfeID = spGridFct->local_finite_element_id(m_fct);
418 0 : };
419 :
420 0 : virtual bool continuous() const
421 : {
422 0 : return false;
423 : }
424 :
425 : template <int refDim>
426 0 : void eval_and_deriv(MathVector<dim> vValue[],
427 : const MathVector<dim> vGlobIP[],
428 : number time, int si,
429 : GridObject* elem,
430 : const MathVector<dim> vCornerCoords[],
431 : const MathVector<refDim> vLocIP[],
432 : const size_t nip,
433 : LocalVector* u,
434 : bool bDeriv,
435 : int s,
436 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[],
437 : const MathMatrix<refDim, dim>* vJT = NULL) const
438 : {
439 : // reference object id
440 0 : const ReferenceObjectID roid = elem->reference_object_id();
441 :
442 : // get reference element mapping by reference object id
443 0 : std::vector<MathMatrix<refDim, dim> > vJTTmp(nip);
444 0 : if(vJT == NULL){
445 : try{
446 : DimReferenceMapping<refDim, dim>& mapping
447 : = ReferenceMappingProvider::get<refDim, dim>(roid, vCornerCoords);
448 :
449 : // compute transformation matrices
450 0 : mapping.jacobian_transposed(&(vJTTmp[0]), vLocIP, nip);
451 :
452 : // store tmp Gradient
453 : vJT = &(vJTTmp[0]);
454 0 : }UG_CATCH_THROW("GridFunctionGradientData: failed.");
455 : }
456 :
457 : // get trial space
458 : try{
459 : const LocalShapeFunctionSet<refDim>& rTrialSpace =
460 0 : LocalFiniteElementProvider::get<refDim>(roid, m_lfeID);
461 :
462 : // storage for shape function at ip
463 : std::vector<MathVector<refDim> > vLocGrad;
464 : MathVector<refDim> locGrad;
465 :
466 : // Reference Mapping
467 : MathMatrix<dim, refDim> JTInv;
468 :
469 : // loop ips
470 0 : for(size_t ip = 0; ip < nip; ++ip)
471 : {
472 : // evaluate at shapes at ip
473 0 : rTrialSpace.grads(vLocGrad, vLocIP[ip]);
474 :
475 : // get multiindices of element
476 : std::vector<DoFIndex > ind;
477 0 : m_spGridFct->dof_indices(elem, m_fct, ind);
478 :
479 : // compute grad at ip
480 : VecSet(locGrad, 0.0);
481 0 : for(size_t sh = 0; sh < vLocGrad.size(); ++sh)
482 : {
483 0 : const number valSH = DoFRef( *m_spGridFct, ind[sh]);
484 : VecScaleAppend(locGrad, valSH, vLocGrad[sh]);
485 : }
486 :
487 0 : RightInverse (JTInv, vJT[ip]);
488 0 : MatVecMult(vValue[ip], JTInv, locGrad);
489 : }
490 0 : }
491 0 : UG_CATCH_THROW("GridFunctionNumberData: Shape Function Set missing for"
492 : " Reference Object: "<<roid<<", Trial Space: "
493 : <<m_lfeID<<", refDim="<<refDim);
494 :
495 0 : if(bDeriv)
496 0 : UG_THROW("Not implemented.");
497 0 : }
498 : };
499 :
500 :
501 : /**
502 : * \brief Retrieve component of gradient of GridFunction
503 : * \details Helper construct to retrieve specific vector element of the gradient
504 : * of a GridFunction inside loops (e.g. integrals) over that GridFunction.
505 : * \code{.lua}
506 : -- integrate the second component of the gradient of the GridFunction u
507 : Integrate( GridFunctionGradientComponentData(u, "c", 2), u )
508 : * \endcode
509 : */
510 : template <typename TGridFunction>
511 : class GridFunctionGradientComponentData
512 : : public StdDependentUserData<GridFunctionGradientComponentData<TGridFunction>,
513 : number, TGridFunction::dim>
514 : {
515 : public:
516 : /// World dimension of GridFunction m_spGridFct
517 : static const int dim = TGridFunction::dim;
518 :
519 : private:
520 : /// GridFunction to loop over
521 : SmartPtr<TGridFunction> m_spGridFct;
522 :
523 : /// Component ID of function to be used
524 : size_t m_fct;
525 :
526 : /// Component index of gradient to return (0-based)
527 : size_t m_component;
528 :
529 : /// Local Finite Element ID
530 : LFEID m_lfeID;
531 :
532 : public:
533 : /**
534 : * \brief Constructor
535 : * \param[in] spGridFct GridFunction to loop over
536 : * \param[in] cmp Name of the GridFunction's function to calculate gradient of
537 : * \param[in] component Index of gradient vector component to return (1-based)
538 : */
539 0 : GridFunctionGradientComponentData( SmartPtr<TGridFunction> spGridFct,
540 : const char* cmp,
541 : size_t component /* 1-based */ )
542 0 : : m_spGridFct( spGridFct )
543 : {
544 0 : this->set_functions(cmp);
545 :
546 : // check validity of component index
547 0 : if ( component > static_cast<size_t>(dim) && component > 0 ) {
548 0 : UG_THROW( "GridFunctionGradientComponentData: Requested component index "
549 : << component << " out of bounds [1," << dim << "]." );
550 : }
551 0 : m_component = component - 1;
552 :
553 : // get function id of name
554 0 : m_fct = spGridFct->fct_id_by_name( cmp );
555 :
556 : // check that function exists
557 0 : if( m_fct >= spGridFct->num_fct() ) {
558 0 : UG_THROW( "GridFunctionGradientComponentData: Function space does not contain"
559 : " a function with name " << cmp << "." );
560 : }
561 :
562 : // local finite element id
563 0 : m_lfeID = spGridFct->local_finite_element_id( m_fct );
564 0 : };
565 :
566 0 : virtual bool continuous() const
567 : {
568 0 : return false;
569 : }
570 :
571 : /**
572 : * \param[out] vValue Array of the <tt>nip</tt> gradient components
573 : */
574 : template <int refDim>
575 0 : void eval_and_deriv(number vValue[],
576 : const MathVector<dim> vGlobIP[],
577 : number time, int si,
578 : GridObject* elem,
579 : const MathVector<dim> vCornerCoords[],
580 : const MathVector<refDim> vLocIP[],
581 : const size_t nip,
582 : LocalVector* u,
583 : bool bDeriv,
584 : int s,
585 : std::vector<std::vector<number> > vvvDeriv[],
586 : const MathMatrix<refDim, dim>* vJT = NULL) const
587 : {
588 : // reference object id
589 0 : const ReferenceObjectID roid = elem->reference_object_id();
590 :
591 : // get reference element mapping by reference object id
592 0 : std::vector<MathMatrix<refDim, dim> > vJTTmp( nip );
593 0 : if( vJT == NULL ) {
594 : try{
595 : DimReferenceMapping<refDim, dim>& mapping
596 : = ReferenceMappingProvider::get< refDim, dim >( roid, vCornerCoords );
597 :
598 : // compute transformation matrices
599 0 : mapping.jacobian_transposed( &(vJTTmp[0]), vLocIP, nip );
600 :
601 : // store tmp Gradient
602 : vJT = &(vJTTmp[0]);
603 0 : } UG_CATCH_THROW( "GridFunctionGradientComponentData: failed.");
604 : }
605 :
606 : // get trial space
607 : try {
608 : const LocalShapeFunctionSet<refDim>& rTrialSpace =
609 0 : LocalFiniteElementProvider::get<refDim>( roid, m_lfeID );
610 :
611 : // storage for shape function at ip
612 : std::vector<MathVector<refDim> > vLocGrad;
613 : MathVector<refDim> locGrad;
614 : std::vector<MathVector<dim> > vValueVec;
615 0 : vValueVec.resize( nip );
616 :
617 : // Reference Mapping
618 : MathMatrix<dim, refDim> JTInv;
619 :
620 : // loop ips
621 0 : for( size_t ip = 0; ip < nip; ++ip ) {
622 : // evaluate at shapes at ip
623 0 : rTrialSpace.grads( vLocGrad, vLocIP[ip] );
624 :
625 : // get multiindices of element
626 : std::vector<DoFIndex> ind;
627 0 : m_spGridFct->dof_indices( elem, m_fct, ind );
628 :
629 : // compute grad at ip
630 : VecSet( locGrad, 0.0 );
631 0 : for( size_t sh = 0; sh < vLocGrad.size(); ++sh ) {
632 0 : const number valSH = DoFRef( *m_spGridFct, ind[sh] );
633 : VecScaleAppend( locGrad, valSH, vLocGrad[sh] );
634 : }
635 :
636 0 : RightInverse( JTInv, vJT[ip] );
637 : MatVecMult( vValueVec[ip], JTInv, locGrad );
638 :
639 0 : vValue[ip] = vValueVec[ip][m_component];
640 :
641 0 : if(bDeriv){
642 0 : for( size_t ip = 0; ip < nip; ++ip ) {
643 : UG_ASSERT(vvvDeriv[ip].size() == 1,
644 : "Single component expected, but "<<vvvDeriv[ip].size())
645 : UG_ASSERT(vvvDeriv[ip][0].size() == vLocGrad.size(),
646 : "Wrong number sh: "<<vvvDeriv[ip][0].size()<<", but expected: "<<vLocGrad.size())
647 :
648 0 : for( size_t sh = 0; sh < vLocGrad.size(); ++sh ) {
649 : MatVecMult( vValueVec[ip], JTInv, vLocGrad[sh] );
650 0 : vvvDeriv[ip][0][sh] = vValueVec[ip][m_component];
651 : }
652 : }
653 : }
654 : }
655 0 : }
656 0 : UG_CATCH_THROW("GridFunctionNumberData: Shape Function Set missing for"
657 : " Reference Object: "<<roid<<", Trial Space: "
658 : <<m_lfeID<<", refDim="<<refDim);
659 0 : }
660 : };
661 :
662 : } // end namespace ug
663 :
664 : #endif /* __H__UG__LIB_DISC__FUNCTION_SPACE__GRID_FUNCTION_USER_DATA__ */
|
