Line data Source code
1 : /*
2 : * Copyright (c) 2013-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 "lagrange_local_dof.h"
34 : #include "common/util/provider.h"
35 : #include "lib_disc/common/multi_index.h"
36 :
37 : namespace ug{
38 :
39 : ///////////////////////////////////////////////////////////////////////////////
40 : // Help Functions to create LocalDoFs
41 : ///////////////////////////////////////////////////////////////////////////////
42 :
43 19 : void SetLagrangeVertexLocalDoFs(std::vector<LocalDoF>& vLocalDoF,
44 : const ReferenceElement& rRef,
45 : const size_t p)
46 : {
47 : // loop all vertices
48 105 : for(size_t i = 0; i< rRef.num(0); ++i)
49 86 : vLocalDoF.push_back(LocalDoF(0, i, 0));
50 19 : }
51 :
52 19 : void SetLagrangeEdgeLocalDoFs(std::vector<LocalDoF>& vLocalDoF,
53 : const ReferenceElement& rRef,
54 : const size_t p)
55 : {
56 : // only for 2d,3d elems we do something
57 19 : if(rRef.dimension() < 1) return;
58 :
59 : // loop all edges
60 132 : for(size_t e = 0; e< rRef.num(1); ++e)
61 : {
62 : // add dofs on the edge
63 218 : for(size_t i = 1; i < p; ++i)
64 : {
65 : // set: dim=1, id=e, offset=i-1
66 105 : vLocalDoF.push_back(LocalDoF(1, e, i-1));
67 : }
68 : }
69 : }
70 :
71 19 : void SetLagrangeFaceLocalDoFs(std::vector<LocalDoF>& vLocalDoF,
72 : const ReferenceElement& rRef,
73 : const size_t p)
74 : {
75 : // only for 2d,3d elems we do something
76 19 : if(rRef.dimension() < 2) return;
77 :
78 : // add dof on quadrilateral
79 72 : for(size_t f = 0; f< rRef.num(2); ++f)
80 : {
81 : // reset counter
82 : size_t cnt = 0;
83 :
84 : // loop 'y'-direction
85 107 : for(size_t j = 1; j < p; ++j)
86 : {
87 : // for a quadrilateral we have a quadratic loop, but for a
88 : // triangle we need to stop at the diagonal
89 51 : const size_t off = ((rRef.num(2,f,0)==3) ? j : 0);
90 :
91 : // loop 'x'-direction
92 108 : for(size_t i = 1; i < p-off; ++i)
93 : {
94 : // set: dim=2, id=f, offset=cnt
95 57 : vLocalDoF.push_back(LocalDoF(2, f, cnt++));
96 : }
97 : }
98 : }
99 : }
100 :
101 19 : void SetLagrangeVolumeLocalDoFs(std::vector<LocalDoF>& vLocalDoF,
102 : const ReferenceElement& rRef,
103 : const size_t p)
104 : {
105 : // only for 3d elems we do something
106 19 : if(rRef.dimension() < 3) return;
107 :
108 : // get type of reference element
109 : const ReferenceObjectID roid = rRef.roid();
110 :
111 : // handle elems
112 : size_t cnt = 0;
113 10 : switch(roid)
114 : {
115 : case ROID_TETRAHEDRON:
116 6 : for(size_t m2 = 1; m2 < p; ++m2)
117 4 : for(size_t m1 = 1; m1 < p-m2; ++m1)
118 1 : for(size_t m0 = 1; m0 < p-m2-m1; ++m0)
119 : {
120 : // set: dim=2, id=0, offset=i
121 0 : vLocalDoF.push_back(LocalDoF(3, 0, cnt++));
122 : }
123 : break;
124 :
125 : case ROID_PYRAMID:
126 : //\todo:order dofs
127 : {
128 : size_t numInnerDoF = 0;
129 1 : for(int i=1; i <= (int)p -2; ++i) numInnerDoF += i*i;
130 :
131 1 : for(size_t i = 0; i < numInnerDoF; ++i)
132 0 : vLocalDoF.push_back(LocalDoF(3, 0, i));
133 : }
134 : break;
135 :
136 : case ROID_PRISM:
137 6 : for(size_t m2 = 1; m2 < p; ++m2)
138 8 : for(size_t m1 = 1; m1 < p; ++m1)
139 7 : for(size_t m0 = 1; m0 < p-m1; ++m0)
140 : {
141 : // set: dim=2, id=0, offset=i
142 2 : vLocalDoF.push_back(LocalDoF(3, 0, cnt++));
143 : }
144 : break;
145 :
146 : case ROID_HEXAHEDRON:
147 6 : for(size_t m2 = 1; m2 < p; ++m2)
148 8 : for(size_t m1 = 1; m1 < p; ++m1)
149 14 : for(size_t m0 = 1; m0 < p; ++m0)
150 : {
151 : // set: dim=2, id=0, offset=i
152 9 : vLocalDoF.push_back(LocalDoF(3, 0, cnt++));
153 : }
154 : break;
155 :
156 0 : case ROID_OCTAHEDRON:
157 : {
158 0 : if(p != 1)
159 0 : UG_THROW("SetLagrangeVolumeLocalDoFs: Octahedral elements only implemented for order p = 1.");
160 : }
161 : break;
162 :
163 0 : default: UG_THROW("SetLagrangeVolumeLocalDoFs: Missing 3d mapping "
164 : "for type '"<<roid<<"'.");
165 : }
166 : }
167 :
168 19 : void SetLagrangeLocalDoFs( std::vector<LocalDoF>& vLocalDoF,
169 : const ReferenceElement& rRef,
170 : const size_t p)
171 : {
172 19 : SetLagrangeVertexLocalDoFs(vLocalDoF, rRef, p);
173 19 : SetLagrangeEdgeLocalDoFs(vLocalDoF, rRef, p);
174 19 : SetLagrangeFaceLocalDoFs(vLocalDoF, rRef, p);
175 19 : SetLagrangeVolumeLocalDoFs(vLocalDoF, rRef, p);
176 :
177 19 : if(vLocalDoF.size() != LagrangeNumDoFs(rRef.roid(), p))
178 0 : UG_THROW("Wrong number of LocalDoFs ("<<vLocalDoF.size()<<") distributed, "
179 : "correct is "<<LagrangeNumDoFs(rRef.roid(), p));
180 19 : }
181 :
182 1 : size_t GetNumberOfDoFsOfPyramid(int p)
183 : {
184 1 : if(p <= 0) return 0;
185 : if(p == 0) return 1;
186 1 : if(p == 1) return 5;
187 0 : else return GetNumberOfDoFsOfPyramid(p-1) + (p+1)*(p+1);
188 : }
189 :
190 0 : size_t LagrangeNumDoFOnSub(const ReferenceObjectID elem,
191 : const ReferenceObjectID sub, const size_t p)
192 : {
193 0 : switch(elem){
194 0 : case ROID_VERTEX:
195 0 : if(sub == ROID_VERTEX) return 1;
196 0 : else return 0;
197 0 : case ROID_EDGE:
198 0 : if(sub == ROID_VERTEX) return 1;
199 0 : else if(sub == ROID_EDGE) return p-1;
200 : else return 0;
201 0 : case ROID_TRIANGLE:
202 0 : if(sub == ROID_VERTEX) return 1;
203 0 : if(sub == ROID_EDGE) return (p-1);
204 0 : if(sub == ROID_TRIANGLE) return ((p>2) ? BinomCoeff(p-1, p-3) : 0);
205 : else return 0;
206 0 : case ROID_QUADRILATERAL:
207 0 : if(sub == ROID_VERTEX) return 1;
208 0 : if(sub == ROID_EDGE) return (p-1);
209 0 : if(sub == ROID_QUADRILATERAL) return (p-1)*(p-1);
210 : else return 0;
211 0 : case ROID_TETRAHEDRON:
212 0 : if(sub == ROID_VERTEX) return 1;
213 0 : if(sub == ROID_EDGE) return (p-1);
214 : // same as for a 2d triangle of order p-3
215 0 : if(sub == ROID_TRIANGLE) return ((p>2) ? BinomCoeff(p-1, p-3) : 0);
216 : // same as for a 3d tetrahedron of order p-4
217 0 : if(sub == ROID_TETRAHEDRON) return ((p>3) ? BinomCoeff(p-1, p-4) : 0);
218 : else return 0;
219 0 : case ROID_PRISM:
220 : if(sub == ROID_VERTEX) return 1;
221 0 : if(sub == ROID_EDGE) return (p-1);
222 : // same as for a 2d triangle of order p-3
223 0 : if(sub == ROID_TRIANGLE) return ((p>2) ? BinomCoeff(p-1, p-3) : 0);
224 : // same as for a 2d quadrilateral of order p-2
225 0 : if(sub == ROID_QUADRILATERAL) return (p-1)*(p-1);
226 : // same as for a 3d prism of order p-2
227 0 : if(sub == ROID_PRISM) return ((p>2) ? BinomCoeff(p-1, p-3)*(p-1) : 0);
228 0 : else return 0;
229 0 : case ROID_PYRAMID:
230 : if(sub == ROID_VERTEX) return 1;
231 0 : if(sub == ROID_EDGE) return (p-1);
232 : // same as for a 2d triangle of order p-3
233 0 : if(sub == ROID_TRIANGLE) return ((p>2) ? BinomCoeff(p-1, p-3) : 0);
234 : // same as for a 2d quadrilateral of order p-2
235 0 : if(sub == ROID_QUADRILATERAL) return (p-1)*(p-1);
236 : // same as for a 3d pyramid of order p-2
237 0 : if(sub == ROID_PYRAMID) return ((p>2) ? GetNumberOfDoFsOfPyramid(p-3) : 0);
238 0 : else return 0;
239 0 : case ROID_HEXAHEDRON:
240 0 : if(sub == ROID_VERTEX) return 1;
241 0 : if(sub == ROID_EDGE) return (p-1);
242 0 : if(sub == ROID_QUADRILATERAL) return (p-1)*(p-1);
243 0 : if(sub == ROID_HEXAHEDRON) return (p-1)*(p-1)*(p-1);
244 : else return 0;
245 0 : case ROID_OCTAHEDRON:
246 0 : if(p != 1)
247 : {
248 0 : UG_THROW("LagrangeNumDoFOnSub: Octahedral elements only implemented for order p = 1.");
249 : }
250 0 : if(sub == ROID_VERTEX) return 1;
251 0 : else return 0;
252 0 : default: UG_THROW("LagrangeLDS: Invalid ReferenceObjectID: "<<elem);
253 : }
254 : }
255 :
256 19 : size_t LagrangeNumDoFs(const ReferenceObjectID elem, const size_t p)
257 : {
258 19 : switch(elem){
259 : case ROID_VERTEX: return 1;
260 3 : case ROID_EDGE: return p+1;
261 3 : case ROID_TRIANGLE: return BinomCoeff(2 + p, p);
262 3 : case ROID_QUADRILATERAL: return (p+1)*(p+1);
263 3 : case ROID_TETRAHEDRON: return BinomCoeff(3 + p, p);
264 3 : case ROID_PRISM: return BinomCoeff(2+p,p) * (p+1);
265 1 : case ROID_PYRAMID: return GetNumberOfDoFsOfPyramid(p);
266 3 : case ROID_HEXAHEDRON: return (p+1)*(p+1)*(p+1);
267 0 : case ROID_OCTAHEDRON:
268 0 : if(p != 1)
269 : {
270 0 : UG_THROW("LagrangeNumDoFs: Octahedral elements only implemented for order p = 1.");
271 : }
272 : else
273 : return 6;
274 0 : default: UG_THROW("LagrangeLDS: Invalid ReferenceObjectID: "<<elem);
275 : }
276 : }
277 :
278 :
279 : ///////////////////////////////////////////////////////////////////////////////
280 : // LagrangeLDS
281 : ///////////////////////////////////////////////////////////////////////////////
282 :
283 : template <typename TRefElem>
284 19 : LagrangeLDS<TRefElem>::LagrangeLDS(size_t order)
285 : {
286 19 : set_order(order);
287 19 : }
288 :
289 : template <typename TRefElem>
290 19 : void LagrangeLDS<TRefElem>::set_order(size_t order)
291 : {
292 19 : p = order;
293 : m_vLocalDoF.clear();
294 19 : SetLagrangeLocalDoFs(m_vLocalDoF, Provider<TRefElem>::get(), p);
295 19 : }
296 :
297 : template <typename TRefElem>
298 0 : size_t LagrangeLDS<TRefElem>::num_dof(ReferenceObjectID type) const
299 : {
300 0 : return LagrangeNumDoFOnSub(roid(), type, p);
301 : }
302 :
303 : template class LagrangeLDS<ReferenceVertex>;
304 : template class LagrangeLDS<ReferenceEdge>;
305 : template class LagrangeLDS<ReferenceTriangle>;
306 : template class LagrangeLDS<ReferenceQuadrilateral>;
307 : template class LagrangeLDS<ReferenceTetrahedron>;
308 : template class LagrangeLDS<ReferencePrism>;
309 : template class LagrangeLDS<ReferencePyramid>;
310 : template class LagrangeLDS<ReferenceHexahedron>;
311 : template class LagrangeLDS<ReferenceOctahedron>;
312 :
313 : } // end namespace ug
|
