Line data Source code
1 : /*
2 : * Copyright (c) 2010-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 : #include "local_finite_element_provider.h"
34 : #include "lib_disc/reference_element/reference_element_util.h"
35 : #include "lib_disc/reference_element/reference_mapping_provider.h"
36 :
37 : // include spaces
38 : #include "lagrange/lagrangep1.h"
39 : #include "lagrange/lagrange.h"
40 : #include "crouzeix-raviart/crouzeix_raviart.h"
41 : #include "piecewise_constant/piecewise_constant.h"
42 : #include "mini/mini.h"
43 : #include "nedelec/nedelec.h"
44 :
45 :
46 : namespace ug{
47 :
48 : DebugID DID_LOCAL_FINITE_ELEMENT_PROVIDER("LOCAL_FINITE_ELEMENT_PROVIDER");
49 :
50 : ////////////////////////////////////////////////////////////////////////////////
51 : // Wrapper for LocalShapeFunctionSets
52 : ////////////////////////////////////////////////////////////////////////////////
53 :
54 : /// wrapper class implementing the LocalShapeFunctionSet interface
55 : /**
56 : * This class wrappes a class passed by the template argument into the
57 : * virtual ILocalShapeFunctionSet interface and makes it thus usable in that
58 : * context on the price of virtual functions.
59 : *
60 : * \tparam TImpl Implementation of a Local Shape Function Set
61 : */
62 : template <typename TImpl>
63 : class LocalShapeFunctionSetWrapper
64 : : public LocalShapeFunctionSet<TImpl::dim,
65 : typename TImpl::shape_type,
66 : typename TImpl::grad_type>,
67 : public TImpl
68 : {
69 : /// Implementation
70 : typedef TImpl ImplType;
71 :
72 : public:
73 : /// reference element dimension
74 : static const int dim = TImpl::dim;
75 :
76 : /// Shape type
77 : typedef typename ImplType::shape_type shape_type;
78 :
79 : /// Gradient type
80 : typedef typename ImplType::grad_type grad_type;
81 :
82 : public:
83 : /// constructor
84 0 : LocalShapeFunctionSetWrapper(){}
85 :
86 : public:
87 : /// \copydoc ug::LocalDoFSet::type()
88 0 : virtual ReferenceObjectID roid() const {return ImplType::roid();}
89 :
90 : /// \copydoc ug::LocalDoFSet::num_sh()
91 0 : virtual size_t num_sh() const {return ImplType::num_sh();}
92 :
93 : /// \copydoc ug::LocalDoFSet::num_dof()
94 0 : virtual size_t num_dof(ReferenceObjectID roid) const {return ImplType::num_dof(roid);}
95 :
96 : /// \copydoc ug::LocalDoFSet::local_dof()
97 0 : virtual const LocalDoF& local_dof(size_t dof) const {return ImplType::local_dof(dof);}
98 :
99 : /// \copydoc ug::LocalDoFSet::exact_position_available()
100 0 : virtual bool exact_position_available() const {return ImplType::exact_position_available();}
101 :
102 : /// \copydoc ug::LocalShapeFunctionSet::position()
103 0 : virtual bool position(size_t i, MathVector<dim>& pos) const
104 : {
105 0 : return ImplType::position(i, pos);
106 : }
107 :
108 : public:
109 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
110 0 : virtual bool continuous() const {return ImplType::continuous();}
111 :
112 : /// \copydoc ug::LocalShapeFunctionSet::shape()
113 0 : virtual shape_type shape(size_t i, const MathVector<dim>& x) const
114 : {
115 0 : return ImplType::shape(i, x);
116 : }
117 :
118 : /// \copydoc ug::LocalShapeFunctionSet::shape()
119 0 : virtual void shape(shape_type& s, size_t i, const MathVector<dim>& x) const
120 : {
121 : typedef BaseLSFS<TImpl, TImpl::dim,
122 : typename TImpl::shape_type,
123 : typename TImpl::grad_type> baseType;
124 :
125 0 : baseType::shape(s, i, x);
126 0 : }
127 :
128 : /// \copydoc ug::LocalShapeFunctionSet::shapes()
129 0 : virtual void shapes(shape_type* vShape, const MathVector<dim>& x) const
130 : {
131 0 : ImplType::shapes(vShape, x);
132 0 : }
133 :
134 : /// \copydoc ug::LocalShapeFunctionSet::shapes()
135 0 : virtual void shapes(std::vector<shape_type>& vShape, const MathVector<dim>& x) const
136 : {
137 0 : ImplType::shapes(vShape, x);
138 0 : }
139 :
140 : /// \copydoc ug::LocalShapeFunctionSet::shapes()
141 0 : virtual void shapes(std::vector<std::vector<shape_type> >& vvShape,
142 : const std::vector<MathVector<dim> >& vLocPos) const
143 : {
144 0 : ImplType::shapes(vvShape, vLocPos);
145 0 : }
146 :
147 : /// \copydoc ug::LocalShapeFunctionSet::grad()
148 0 : virtual void grad(grad_type& g, size_t i, const MathVector<dim>& x) const
149 : {
150 0 : ImplType::grad(g, i, x);
151 0 : }
152 :
153 : /// \copydoc ug::LocalShapeFunctionSet::grads()
154 0 : virtual void grads(grad_type* vGrad, const MathVector<dim>& x) const
155 : {
156 0 : ImplType::grads(vGrad, x);
157 0 : }
158 :
159 : /// \copydoc ug::LocalShapeFunctionSet::grads()
160 0 : virtual void grads(std::vector<grad_type>& vGrad, const MathVector<dim>& x) const
161 : {
162 0 : ImplType::grads(vGrad, x);
163 0 : }
164 :
165 : /// \copydoc ug::LocalShapeFunctionSet::grads()
166 0 : virtual void grads(std::vector<std::vector<grad_type> >& vvGrad,
167 : const std::vector<MathVector<dim> >& vLocPos) const
168 : {
169 0 : ImplType::grads(vvGrad, vLocPos);
170 0 : }
171 : };
172 :
173 : ////////////////////////////////////////////////////////////////////////////////
174 : // SubLocalDoFSet
175 : ////////////////////////////////////////////////////////////////////////////////
176 :
177 : /**
178 : * Intersection of local dofs
179 : */
180 : template <int TDim>
181 : class SubLocalDoFSet : public DimLocalDoFSet<TDim>
182 : {
183 : static const int dim = TDim;
184 :
185 : public:
186 : template <int setDim>
187 0 : SubLocalDoFSet(const ReferenceObjectID roid,
188 : const DimLocalDoFSet<setDim>& set)
189 0 : : m_roid(roid),
190 0 : m_bExactPos(set.exact_position_available()),
191 0 : m_bInit(false),
192 0 : m_vNumDoF(NUM_REFERENCE_OBJECTS, 0)
193 : {
194 0 : if(ReferenceElementDimension(roid) != dim)
195 0 : UG_THROW("SubLocalDoFSet: templated Reference Dimension "<<dim<<
196 : " does not match Reference Element dimension of "<<roid);
197 :
198 : const DimReferenceElement<setDim>& setRefElem =
199 0 : ReferenceElementProvider::get<setDim>(set.roid());
200 :
201 : // loop all subs of roid type
202 0 : for(size_t id = 0; id < setRefElem.num(dim); ++id){
203 0 : if(setRefElem.roid(dim, id) != m_roid) continue;
204 :
205 : // get mapping to sub
206 : std::vector<MathVector<setDim> > vCorner;
207 0 : for(size_t co = 0; co < setRefElem.num(dim, id, 0); ++co){
208 0 : vCorner.push_back(setRefElem.corner(setRefElem.id(dim, id, 0, co)));
209 : }
210 : const DimReferenceMapping<dim, setDim>& map =
211 0 : ReferenceMappingProvider::get<dim, setDim>(m_roid, vCorner);
212 :
213 0 : std::vector<size_t> vNumDoF(NUM_REFERENCE_OBJECTS, 0);
214 : std::vector<LocalDoF> vLocalDoF;
215 : std::vector<MathVector<dim> > vLocalPos;
216 :
217 : // loop subselements of sub
218 0 : for(int subDim = 0; subDim <= dim; ++subDim){
219 :
220 0 : for(size_t i = 0; i < setRefElem.num(dim, id, subDim); ++i){
221 0 : const size_t subID = setRefElem.id(dim, id, subDim, i);
222 :
223 : const ReferenceObjectID sroid = setRefElem.roid(subDim, subID);
224 0 : vNumDoF[sroid] = set.num_dof(sroid);
225 :
226 : // get local DoFs on sub
227 : std::vector<size_t> vLocalDoFID;
228 0 : for(size_t dof = 0; dof < set.num_dof(); ++dof){
229 0 : const LocalDoF& localDoF = set.local_dof(dof);
230 0 : if(localDoF.dim() != subDim) continue;
231 0 : if(localDoF.id() != subID) continue;
232 :
233 0 : vLocalDoFID.push_back(dof);
234 : }
235 :
236 : // loop sub DoFs in correct order
237 0 : for(size_t offset = 0; offset < vLocalDoFID.size(); ++offset){
238 :
239 : const LocalDoF* pLocalDoF = NULL;
240 : size_t localDoFID = (size_t)-1;
241 0 : for(size_t dof = 0; dof < vLocalDoFID.size(); ++dof){
242 0 : if(set.local_dof(vLocalDoFID[dof]).offset() == offset){
243 0 : pLocalDoF = &set.local_dof(vLocalDoFID[dof]);
244 0 : localDoFID = vLocalDoFID[dof];
245 : }
246 : }
247 :
248 0 : if(pLocalDoF == NULL)
249 0 : UG_THROW("SubLocalDoFSet: Cannot find local dof "
250 : "with offset "<<offset<<", but must "
251 : "exist. Check Implementation of LocalDoFSet ");
252 :
253 : // map local dof
254 : MathVector<dim> locPos(0.0);
255 : MathVector<setDim> globPos;
256 0 : set.position(localDoFID, globPos);
257 :
258 : try{
259 0 : map.global_to_local(locPos, globPos);
260 : }
261 0 : catch(UGError& err){
262 0 : std::stringstream ss;
263 : ss<<"SubLocalDoFSet: Cannot find local position for "
264 0 : "global Position "<<globPos<<" of the "<<localDoFID<<
265 0 : "'th LocalDoF on "<<set.roid()<<
266 0 : " for the "<<id<<"'th "<<roid<<" and its "
267 0 : "Sub-"<<sroid<<" with id "<<subID<<" with "
268 : "corners:\n";
269 0 : for(size_t i = 0; i < vCorner.size(); ++i)
270 0 : ss << " "<<i<<": "<<vCorner[i] << "\n";
271 :
272 0 : err.push_msg(ss.str(),__FILE__,__LINE__);
273 0 : throw(err);
274 0 : }
275 :
276 : // add
277 0 : vLocalDoF.push_back(LocalDoF(subDim, i, offset));
278 0 : vLocalPos.push_back(locPos);
279 : }
280 : }
281 : }
282 :
283 0 : this->set(vNumDoF, vLocalDoF, vLocalPos);
284 : }
285 0 : }
286 :
287 : public:
288 0 : virtual ReferenceObjectID roid() const {return m_roid;};
289 0 : virtual size_t num_dof(ReferenceObjectID roid) const {return m_vNumDoF[roid];}
290 0 : virtual size_t num_sh() const {return m_vLocalDoF.size();}
291 0 : virtual const LocalDoF& local_dof(size_t dof) const {return m_vLocalDoF[dof];}
292 0 : virtual bool exact_position_available() const {return m_bExactPos;}
293 0 : virtual bool position(size_t i, MathVector<TDim>& pos) const
294 : {
295 : pos = m_vLocalPos[i];
296 0 : return exact_position_available();
297 : }
298 :
299 : protected:
300 0 : void set(const std::vector<size_t>& vNumDoF,
301 : const std::vector<LocalDoF>& vLocalDoF,
302 : const std::vector<MathVector<TDim> >& vLocalPos)
303 : {
304 0 : if(m_bInit){
305 0 : if(vNumDoF != m_vNumDoF) UG_THROW("NumDoF mismatch");
306 0 : if(vLocalDoF != m_vLocalDoF) UG_THROW("vLocalDoF mismatch");
307 0 : if(vLocalPos != m_vLocalPos) UG_THROW("vLocalPos mismatch");
308 : }
309 : else{
310 0 : m_vNumDoF = vNumDoF;
311 0 : m_vLocalDoF = vLocalDoF;
312 0 : m_vLocalPos = vLocalPos;
313 : }
314 0 : }
315 :
316 :
317 : protected:
318 : const ReferenceObjectID m_roid; ///< Reference ID this DoF Set is for
319 : const bool m_bExactPos;
320 :
321 : bool m_bInit;
322 : std::vector<size_t> m_vNumDoF;
323 : std::vector<LocalDoF> m_vLocalDoF; ///< Local DoFs of this set
324 : std::vector<MathVector<TDim> > m_vLocalPos; ///< Local Positions of DoFs
325 : };
326 :
327 : ////////////////////////////////////////////////////////////////////////////////
328 : // Provider for LocalShapeFunctionSets
329 : ////////////////////////////////////////////////////////////////////////////////
330 :
331 : template <typename TRefElem>
332 0 : void LocalFiniteElementProvider::create_lagrange_set(const LFEID& id)
333 : {
334 : // reference object id
335 : static const ReferenceObjectID roid = TRefElem::REFERENCE_OBJECT_ID;
336 : static const int dim = TRefElem::dim;
337 :
338 : // only order >= 1 available
339 0 : if(id.order() < 1) return;
340 :
341 : // if refdim == dim create the space
342 0 : if(dim == id.dim()){
343 : ConstSmartPtr<LocalShapeFunctionSet<dim> > set;
344 0 : switch(id.order()){
345 0 : case 1: set = ConstSmartPtr<LocalShapeFunctionSet<dim> >
346 0 : (new LocalShapeFunctionSetWrapper<LagrangeP1<TRefElem> >);
347 0 : break;
348 0 : case 2: set = ConstSmartPtr<LocalShapeFunctionSet<dim> >
349 0 : (new LocalShapeFunctionSetWrapper<LagrangeLSFS<TRefElem,2> >);
350 0 : break;
351 : break;
352 0 : case 3: set = ConstSmartPtr<LocalShapeFunctionSet<dim> >
353 0 : (new LocalShapeFunctionSetWrapper<LagrangeLSFS<TRefElem,3> >);
354 0 : break;
355 : break;
356 0 : case 4: set = ConstSmartPtr<LocalShapeFunctionSet<dim> >
357 0 : (new LocalShapeFunctionSetWrapper<LagrangeLSFS<TRefElem,4> >);
358 0 : break;
359 : break;
360 :
361 0 : default:{
362 : SmartPtr<LocalShapeFunctionSetWrapper<FlexLagrangeLSFS<TRefElem> > >
363 0 : sSetFlexLagrange(new LocalShapeFunctionSetWrapper<FlexLagrangeLSFS<TRefElem> >);
364 0 : sSetFlexLagrange->set_order(id.order());
365 0 : set = sSetFlexLagrange;
366 : }
367 : }
368 0 : register_set(id, set);
369 : return;
370 : }
371 : // if refdim < dim, the restriction of to the subelement is the lagrange space
372 : // of the refdim
373 0 : else if (dim < id.dim()){
374 : // get (and maybe create) lagrange space for dim == refdim (always exist)
375 : ConstSmartPtr<LocalShapeFunctionSet<dim> > set =
376 0 : getptr<dim>(roid, LFEID(id.type(), dim, id.order()));
377 :
378 : // register this space as space on subelement
379 0 : register_set(id, set);
380 : }
381 :
382 : // else: for refdim > dim there exist no restrictions
383 : }
384 :
385 0 : void LocalFiniteElementProvider::
386 : create_lagrange_set(ReferenceObjectID roid, const LFEID& id)
387 : {
388 : UG_DLOG(DID_LOCAL_FINITE_ELEMENT_PROVIDER, 2, ">>OCT_DISC_DEBUG: " << "local_finite_element_provider.cpp: " << "create_lagrange_set(): " << roid << std::endl);
389 0 : switch(roid){
390 0 : case ROID_VERTEX: create_lagrange_set<ReferenceVertex>(id); return;
391 0 : case ROID_EDGE: create_lagrange_set<ReferenceEdge>(id); return;
392 0 : case ROID_TRIANGLE: create_lagrange_set<ReferenceTriangle>(id); return;
393 0 : case ROID_QUADRILATERAL:create_lagrange_set<ReferenceQuadrilateral>(id); return;
394 0 : case ROID_TETRAHEDRON: create_lagrange_set<ReferenceTetrahedron>(id); return;
395 0 : case ROID_PRISM: create_lagrange_set<ReferencePrism>(id); return;
396 0 : case ROID_HEXAHEDRON: create_lagrange_set<ReferenceHexahedron>(id); return;
397 : case ROID_PYRAMID:
398 : // only space available for order 1
399 0 : if(id.order() != 1)
400 : return;
401 : else
402 : {
403 0 : register_set(LFEID(LFEID::LAGRANGE, 3, 1),
404 0 : ConstSmartPtr<LocalShapeFunctionSet<3> >
405 0 : (new LocalShapeFunctionSetWrapper<LagrangeP1<ReferencePyramid> >));
406 : }
407 0 : return;
408 : case ROID_OCTAHEDRON:
409 : // only space available for order 1
410 0 : if(id.order() != 1)
411 : return;
412 : else
413 : {
414 0 : register_set(LFEID(LFEID::LAGRANGE, 3, 1),
415 0 : ConstSmartPtr<LocalShapeFunctionSet<3> >
416 0 : (new LocalShapeFunctionSetWrapper<LagrangeP1<ReferenceOctahedron> >));
417 : }
418 0 : return;
419 : default: return;
420 : }
421 : }
422 :
423 : template <typename TRefElem>
424 0 : void LocalFiniteElementProvider::create_mini_bubble_set(const LFEID& id)
425 : {
426 : // reference object id
427 : static const ReferenceObjectID roid = TRefElem::REFERENCE_OBJECT_ID;
428 : static const int dim = TRefElem::dim;
429 :
430 : // if refdim == dim create the space
431 0 : if(dim == id.dim()){
432 :
433 : if(id == LFEID(LFEID::MINI, dim, 1))
434 0 : register_set(id, ConstSmartPtr<LocalShapeFunctionSet<dim> >(
435 0 : new LocalShapeFunctionSetWrapper<MiniBubbleLSFS<TRefElem> >));
436 0 : return;
437 :
438 : /*
439 : ConstSmartPtr<LocalShapeFunctionSet<dim> > set;
440 :
441 : switch(id.order()){
442 : case 1: set = ConstSmartPtr<LocalShapeFunctionSet<dim> >
443 : (new LocalShapeFunctionSetWrapper<MiniBubbleLSFS<TRefElem> >);
444 : break;
445 : case 2: set = ConstSmartPtr<LocalShapeFunctionSet<dim> >
446 : (new LocalShapeFunctionSetWrapper<MiniBubbleLSFS<TRefElem,2> >);
447 : break;
448 :
449 : }*/
450 :
451 : }
452 : // if refdim < dim, the restriction of to the subelement is the lagrange space
453 : // of the refdim
454 0 : else if (dim < id.dim()){
455 : // get (and maybe create) lagrange space for dim == refdim (always exist)
456 : ConstSmartPtr<LocalShapeFunctionSet<dim> > set =
457 0 : getptr<dim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
458 :
459 : // register this space as space on subelement
460 0 : register_set(id, set);
461 : }
462 :
463 : // else: for refdim > dim there exist no restrictions
464 : }
465 :
466 0 : void LocalFiniteElementProvider::
467 : create_mini_bubble_set(ReferenceObjectID roid, const LFEID& id)
468 : {
469 0 : switch(roid){
470 0 : case ROID_EDGE: create_mini_bubble_set<ReferenceEdge>(id); return;
471 0 : case ROID_TRIANGLE: create_mini_bubble_set<ReferenceTriangle>(id); return;
472 0 : case ROID_QUADRILATERAL:create_mini_bubble_set<ReferenceQuadrilateral>(id); return;
473 0 : case ROID_TETRAHEDRON: create_mini_bubble_set<ReferenceTetrahedron>(id); return;
474 : // no spaces implemented for other 3d elems
475 : default: return;
476 : }
477 : }
478 :
479 : template <typename TRefElem>
480 0 : void LocalFiniteElementProvider::create_nedelec_set(const LFEID& id)
481 : {
482 : static const int dim = TRefElem::dim;
483 : if(id == LFEID(LFEID::NEDELEC, dim, 1))
484 0 : register_set(id, ConstSmartPtr<LocalShapeFunctionSet<dim, MathVector<dim>, MathMatrix<dim,dim> > >(
485 0 : new LocalShapeFunctionSetWrapper<NedelecLSFS<TRefElem> >));
486 0 : return;
487 : }
488 :
489 0 : void LocalFiniteElementProvider::
490 : create_nedelec_set(ReferenceObjectID roid, const LFEID& id)
491 : {
492 0 : switch(roid){
493 0 : case ROID_TRIANGLE: create_nedelec_set<ReferenceTriangle>(id); return;
494 0 : case ROID_TETRAHEDRON: create_nedelec_set<ReferenceTetrahedron>(id); return;
495 : default: return;
496 : }
497 : }
498 :
499 : template <typename TRefElem>
500 0 : void LocalFiniteElementProvider::create_piecewise_constant_set(const LFEID& id)
501 : {
502 : static const int dim = TRefElem::dim;
503 : if(id == LFEID(LFEID::PIECEWISE_CONSTANT, dim, 0))
504 0 : register_set(id, ConstSmartPtr<LocalShapeFunctionSet<dim> >(
505 0 : new LocalShapeFunctionSetWrapper<PiecewiseConstantLSFS<TRefElem> >));
506 0 : return;
507 : }
508 :
509 0 : void LocalFiniteElementProvider::
510 : create_piecewise_constant_set(ReferenceObjectID roid, const LFEID& id)
511 : {
512 0 : switch(roid){
513 0 : case ROID_EDGE: create_piecewise_constant_set<ReferenceEdge>(id); return;
514 0 : case ROID_TRIANGLE: create_piecewise_constant_set<ReferenceTriangle>(id); return;
515 0 : case ROID_QUADRILATERAL:create_piecewise_constant_set<ReferenceQuadrilateral>(id); return;
516 0 : case ROID_TETRAHEDRON: create_piecewise_constant_set<ReferenceTetrahedron>(id); return;
517 0 : case ROID_PRISM: create_piecewise_constant_set<ReferencePrism>(id); return;
518 0 : case ROID_PYRAMID: create_piecewise_constant_set<ReferencePyramid>(id); return;
519 0 : case ROID_HEXAHEDRON: create_piecewise_constant_set<ReferenceHexahedron>(id); return;
520 0 : case ROID_OCTAHEDRON: create_piecewise_constant_set<ReferenceOctahedron>(id); return;
521 : default: return;
522 : }
523 : }
524 :
525 :
526 : template <typename TRefElem>
527 0 : void LocalFiniteElementProvider::create_crouxeiz_raviart_set(const LFEID& id)
528 : {
529 : static const int dim = TRefElem::dim;
530 : if(id == LFEID(LFEID::CROUZEIX_RAVIART, dim, 1))
531 0 : register_set(id, ConstSmartPtr<LocalShapeFunctionSet<dim> >(
532 0 : new LocalShapeFunctionSetWrapper<CrouzeixRaviartLSFS<TRefElem> >));
533 0 : return;
534 : }
535 :
536 0 : void LocalFiniteElementProvider::
537 : create_crouxeiz_raviart_set(ReferenceObjectID roid, const LFEID& id)
538 : {
539 0 : switch(roid){
540 : // no space for dim <= 1
541 0 : case ROID_TRIANGLE: create_crouxeiz_raviart_set<ReferenceTriangle>(id); return;
542 0 : case ROID_QUADRILATERAL:create_crouxeiz_raviart_set<ReferenceQuadrilateral>(id); return;
543 0 : case ROID_TETRAHEDRON: create_crouxeiz_raviart_set<ReferenceTetrahedron>(id); return;
544 0 : case ROID_PRISM: create_crouxeiz_raviart_set<ReferencePrism>(id); return;
545 0 : case ROID_PYRAMID: create_crouxeiz_raviart_set<ReferencePyramid>(id); return;
546 0 : case ROID_HEXAHEDRON: create_crouxeiz_raviart_set<ReferenceHexahedron>(id); return;
547 : default: return;
548 : }
549 : }
550 :
551 0 : void LocalFiniteElementProvider::
552 : create_set(ReferenceObjectID roid, const LFEID& id)
553 : {
554 : try{
555 0 : switch(id.type())
556 : {
557 0 : case LFEID::LAGRANGE: create_lagrange_set(roid, id); break;
558 0 : case LFEID::MINI: create_mini_bubble_set(roid, id); break;
559 0 : case LFEID::NEDELEC: create_nedelec_set(roid, id); break;
560 0 : case LFEID::PIECEWISE_CONSTANT: create_piecewise_constant_set(roid, id); break;
561 0 : case LFEID::CROUZEIX_RAVIART: create_crouxeiz_raviart_set(roid, id); break;
562 : default: return;
563 : }
564 : }
565 0 : UG_CATCH_THROW("LocalFiniteElementProvider: Creation of set "<<id<<
566 : " for "<<roid<<" failed.");
567 : }
568 :
569 0 : void LocalFiniteElementProvider::
570 : create_set(const LFEID& id)
571 : {
572 0 : switch(id.type())
573 : {
574 : // this spaces can create also sub-elem spaces, since they are continuous
575 : case LFEID::LAGRANGE:
576 : case LFEID::MINI:
577 :
578 0 : for(int r = 0; r < NUM_REFERENCE_OBJECTS; ++r){
579 : const ReferenceObjectID roid = (ReferenceObjectID) r;
580 0 : const int dim = ReferenceElementDimension(roid);
581 :
582 0 : if(dim <= id.dim())
583 0 : create_set(roid, id);
584 : }
585 : break;
586 :
587 : // this spaces can not create also sub-elem spaces, since they are discontinuous
588 : case LFEID::PIECEWISE_CONSTANT:
589 : case LFEID::CROUZEIX_RAVIART:
590 : case LFEID::NEDELEC:
591 :
592 0 : for(int r = 0; r < NUM_REFERENCE_OBJECTS; ++r){
593 : const ReferenceObjectID roid = (ReferenceObjectID) r;
594 0 : const int dim = ReferenceElementDimension(roid);
595 :
596 0 : if(dim == id.dim())
597 0 : create_set((ReferenceObjectID)roid, id);
598 : }
599 : break;
600 :
601 : default: return;
602 : }
603 : }
604 :
605 : template <int rdim, int dim>
606 0 : void LocalFiniteElementProvider::
607 : create_sub_dof_set(ReferenceObjectID roid, const LFEID& id)
608 : {
609 0 : if(dim != id.dim())
610 0 : UG_THROW("Dimension must match here, internal error. ("<<dim<<","<<id.dim()<<")");
611 :
612 : // we like to have a DimLocalDoFSet for roid in dimension dim.
613 : // Say roid has dimension rdim. If rdim == dim this must be implemented.
614 : // If rdim < dim, this can be deduced from the implemented patterns for
615 : // elements in dim if and only if some implementation is for a element that
616 : // contains roid as a sub element
617 :
618 0 : for(int r = 0; r < NUM_REFERENCE_OBJECTS; ++r){
619 : const ReferenceObjectID elemRoid = (ReferenceObjectID) r;
620 0 : const int elemDim = ReferenceElementDimension(elemRoid);
621 :
622 : // we only take elements that are in the dimension of the space
623 0 : if(elemDim != dim) continue;
624 :
625 : // try to get the implementation, if not present skip
626 0 : ConstSmartPtr<DimLocalDoFSet<dim> > set = get_dof_ptr<dim>(elemRoid, id);
627 0 : if(set.invalid()) continue;
628 :
629 : // see if element contains the roid
630 : const ReferenceElement& rRefElem = ReferenceElementProvider::get(elemRoid);
631 0 : if(rRefElem.num(roid) == 0) continue;
632 :
633 : try{
634 : // create the sub-dof-pattern
635 0 : ConstSmartPtr<DimLocalDoFSet<rdim> > subSet =
636 0 : ConstSmartPtr<DimLocalDoFSet<rdim> >(new SubLocalDoFSet<rdim>(roid, *set) );
637 :
638 : // try to get an already registerd one
639 0 : ConstSmartPtr<DimLocalDoFSet<rdim> > givenSubSet = get_dof_ptr<rdim>(roid, id, false);
640 :
641 : // if already one set given, check equality
642 0 : if(givenSubSet.valid()){
643 0 : if(*givenSubSet != *subSet)
644 0 : UG_THROW("LocalFiniteElementProvider::create_sub_dof_set:\n"
645 : "Creating DimLocalDoFSet for "<<roid<<" and "<<id<<
646 : ".\nAlready registered Set does not match with computed,"
647 : " but this indicates that the \nLocalDoFSets have not "
648 : "been implemented correctly. Check implementation.\n"
649 : "Sets are: \nGiven:\n"<<*givenSubSet<<"\nvs. New:\n"<<*subSet);
650 : }
651 : else {
652 : // if correct, register set
653 0 : register_set<rdim>(id, subSet);
654 : }
655 : }
656 0 : UG_CATCH_THROW("LocalFiniteElementProvider::create_sub_dof_set: Cannot "
657 : "create SubDoFSet for "<<roid<<" and space "<<id<<" when "
658 : "trying to deduce from LocalDoFSet for "<<elemRoid<<".");
659 : }
660 0 : }
661 :
662 0 : void LocalFiniteElementProvider::
663 : create_dof_set(ReferenceObjectID roid, const LFEID& id)
664 : {
665 : const int dim = id.dim();
666 0 : const int rdim = ReferenceElementDimension(roid);
667 :
668 0 : if(rdim == dim) {
669 0 : create_set(roid, id);
670 0 : return;
671 : }
672 :
673 0 : if(rdim > dim)
674 0 : UG_THROW("Wrong dimension for SubDoFs: "<<rdim<<", "<<dim);
675 :
676 0 : switch(dim){
677 0 : case 1:
678 0 : switch(rdim){
679 0 : case 0: create_sub_dof_set<0, 1>(roid, id); return;
680 : default: return;
681 : }
682 0 : case 2:
683 0 : switch(rdim){
684 0 : case 0: create_sub_dof_set<0, 2>(roid, id); return;
685 0 : case 1: create_sub_dof_set<1, 2>(roid, id); return;
686 : default: return;
687 : }
688 0 : case 3:
689 0 : switch(rdim){
690 0 : case 0: create_sub_dof_set<0, 3>(roid, id); return;
691 0 : case 1: create_sub_dof_set<1, 3>(roid, id); return;
692 0 : case 2: create_sub_dof_set<2, 3>(roid, id); return;
693 : default: return;
694 : }
695 : default: return;
696 : }
697 : }
698 :
699 :
700 0 : LocalFiniteElementProvider::
701 : LocalFiniteElementProvider()
702 0 : {};
703 :
704 0 : LocalFiniteElementProvider::
705 : ~LocalFiniteElementProvider()
706 0 : {};
707 :
708 0 : const LocalDoFSet& LocalFiniteElementProvider::
709 : get_dofs(ReferenceObjectID roid, const LFEID& id, bool bCreate)
710 : {
711 : // init provider and search for identifier
712 : typedef std::map<LFEID, LocalDoFSets> Map;
713 0 : Map::const_iterator iter = inst().m_mLocalDoFSets.find(id);
714 :
715 : // if not found
716 0 : if(iter == m_mLocalDoFSets.end() || (iter->second)[roid].invalid()){
717 0 : if(bCreate){
718 0 : create_dof_set(roid, id);
719 0 : return get_dofs(roid, id, false);
720 : }
721 0 : UG_THROW("LocalFiniteElementProvider: Cannot create LocalDoFSet for finite "
722 : "element type "<<id<<" and "<<roid);
723 : }
724 :
725 : // return dof set
726 : return *((iter->second)[roid]);
727 : }
728 :
729 0 : const CommonLocalDoFSet& LocalFiniteElementProvider::
730 : get_dofs(const LFEID& id, bool bCreate)
731 : {
732 : // init provider and search for identifier
733 : typedef std::map<LFEID, CommonLocalDoFSet> Map;
734 0 : Map::const_iterator iter = inst().m_mCommonDoFSet.find(id);
735 :
736 : // if not found
737 0 : if(iter == m_mCommonDoFSet.end()){
738 0 : if(bCreate) {
739 0 : create_set(id);
740 0 : return get_dofs(id, false);
741 : }
742 0 : UG_THROW("LocalFiniteElementProvider: Cannot create CommonLocalDoFSet for "<<id);
743 : }
744 :
745 : // return the common set
746 0 : return iter->second;
747 : }
748 :
749 0 : bool LocalFiniteElementProvider::continuous(const LFEID& id, bool bCreate)
750 : {
751 : std::map<LFEID, bool>::iterator iter = m_mContSpace.find(id);
752 0 : if(iter == m_mContSpace.end())
753 : {
754 : // try to create the set
755 0 : if(bCreate){
756 0 : create_set(id);
757 0 : return continuous(id, false);
758 : }
759 :
760 0 : UG_THROW("LocalFiniteElementProvider::continuous: no shape function "
761 : "set "<<id<<" registered.");
762 : }
763 :
764 0 : return (*iter).second;
765 : }
766 :
767 0 : void LocalFiniteElementProvider::register_set(const LFEID& id, ConstSmartPtr<LocalDoFSet> set)
768 : {
769 : // reference object id
770 0 : const ReferenceObjectID roid = set->roid();
771 :
772 : // register
773 0 : m_mLocalDoFSets[id][roid] = set;
774 :
775 : // for creation of CommonLocalDoFSet: skip if not the dimension of the space
776 0 : if(set->dim() != id.dim()) return;
777 :
778 : // add this local dof set
779 : try{
780 0 : m_mCommonDoFSet[id].add(*set);
781 : }
782 0 : catch(UGError& err)
783 : {
784 : // write error message
785 0 : std::stringstream ss;
786 : ss<<"LocalFiniteElementProvider::register_set(): "
787 0 : "Cannot build CommonLocalDoFSet for type: "<<id<<" when adding "
788 0 : " Reference element type "<<roid<<".\n"<<
789 0 : "CommonLocalDoFSet is:\n" << m_mCommonDoFSet[id]<<
790 0 : "LocalDoFSet is:\n" << *set;
791 0 : err.push_msg(ss.str(),__FILE__,__LINE__);
792 :
793 : // remove set
794 : m_mCommonDoFSet.erase(id);
795 :
796 0 : throw(err);
797 0 : }
798 : }
799 :
800 : std::map<LFEID, LocalFiniteElementProvider::LocalDoFSets>
801 : LocalFiniteElementProvider::m_mLocalDoFSets =
802 : std::map<LFEID, LocalFiniteElementProvider::LocalDoFSets>();
803 :
804 : std::map<LFEID, CommonLocalDoFSet>
805 : LocalFiniteElementProvider::m_mCommonDoFSet =
806 : std::map<LFEID, CommonLocalDoFSet>();
807 :
808 : std::map<LFEID, bool>
809 : LocalFiniteElementProvider::m_mContSpace =
810 : std::map<LFEID, bool>();
811 :
812 : } // namespace ug
813 :
|
