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 "orientation.h"
34 : #include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h"
35 :
36 : namespace ug{
37 :
38 :
39 : ///////////////////////////////////////////////////////////////////////////////
40 : // Lagrange Offsets
41 : ///////////////////////////////////////////////////////////////////////////////
42 :
43 : /*
44 : * Lagrange DoF Orientation of an Edge:
45 : * If DoFs are assigned to a lower-dimensional edge and we have a degree higher
46 : * than 2 (i.e. more than one DoF on the edge) orientation is required to
47 : * ensure continuity of the shape functions. This means, that each element
48 : * that has the edge as a subelement, must number the dofs on the edge equally
49 : * in global numbering.
50 : *
51 : * The idea is as follows: We induce a global ordering of dofs on the edge by
52 : * using the vertices of the edge itself. We define, that dofs are always
53 : * assigned in a line from the vertex 0 to the vertex 1.
54 : * Now, in the local ordering of dofs on the reference element, the edge may
55 : * have been a different numbering for the corners.
56 : * Thus, we have to distinguish two case:
57 : * a) Orientation matches: we can simply use the usual offset numbering
58 : * b) Orientation mismatches: we have to use the reverse order as offset numbering
59 : */
60 0 : bool OrientationMatches(const EdgeVertices& e1, const EdgeVertices& e2)
61 : {
62 0 : return e1.vertex(0) == e2.vertex(0);
63 : }
64 :
65 0 : void ComputeOrientationOffsetLagrange(std::vector<size_t>& vOrientOffset,
66 : EdgeDescriptor& ed, EdgeVertices* edge, const size_t p)
67 : {
68 0 : vOrientOffset.reserve(p-1);
69 : vOrientOffset.clear();
70 :
71 : // the standard orientation is from co0 -> co1.
72 0 : if(OrientationMatches(ed, *edge))
73 : {
74 0 : for(size_t i = 0; i < p-1; ++i)
75 0 : vOrientOffset.push_back(i);
76 : }
77 : // ... and for reverse order
78 : else
79 : {
80 0 : for(int i = ((int)p) - 2; i >= 0; --i)
81 0 : vOrientOffset.push_back(i);
82 : }
83 0 : }
84 :
85 0 : void ComputeOrientationOffsetLagrange(std::vector<size_t>& vOrientOffset,
86 : Face* face, Edge* edge, size_t nrEdge,
87 : const size_t p)
88 : {
89 0 : EdgeDescriptor ed;
90 0 : face->edge_desc(nrEdge, ed);
91 0 : ComputeOrientationOffsetLagrange(vOrientOffset, ed, edge, p);
92 0 : }
93 :
94 0 : void ComputeOrientationOffsetLagrange(std::vector<size_t>& vOrientOffset,
95 : Volume* vol, Edge* edge, size_t nrEdge,
96 : const size_t p)
97 : {
98 0 : EdgeDescriptor ed;
99 0 : vol->edge_desc(nrEdge, ed);
100 0 : ComputeOrientationOffsetLagrange(vOrientOffset, ed, edge, p);
101 0 : }
102 :
103 : /*
104 : * Lagrange DoF Orientation of a Face:
105 : * If DoFs are assigned to a lower-dimensional face and we have a degree higher
106 : * than 2 (i.e. more than one DoF on the face) orientation is required to
107 : * ensure continuity of the shape functions. This means, that each element
108 : * that has the face as a subelement, must number the dofs on the face equally
109 : * in global numbering.
110 : *
111 : * The idea of the ordering is as follows:
112 : * DoFs are always assigned to the face in a natural order, i.e. in an order such
113 : * that the numbering of the face itself gives the ordering. Now, give a 3d
114 : * element with that face, the reference elements provides us with a numbering
115 : * of the vertices of the face. This numbering must not match the natural ordering
116 : * as present in the face itself.
117 : * Now find the vertex-id (id0) of the face that matches the natural vertex 0 of
118 : * the face itself, and the vertex-id (id1) of the natural vertex 1.
119 : * On the natural face we have a situation like this:
120 : *
121 : * * *--------* ^ j
122 : * | \ | | |
123 : * | \ | | |
124 : * | \ | | |-----> i
125 : * *------ * *--------*
126 : * id0 id1 id0 id1
127 : *
128 : * We define that the DoFs on the face are always numbered in x direction first,
129 : * continuing in the next row in y direction, numbering in x again and
130 : * continuing in y, etc. E.g. this gives (showing only inner dofs):
131 : *
132 : * * *-------* ^ j
133 : * |5 \ | 6 7 8 | |
134 : * |3 4 \ | 3 4 5 | |
135 : * |0 1 2 \ | 0 1 2 | |-----> i
136 : * *-------* *-------*
137 : * id0 id1 id0 id1
138 : * p = 5 p = 4
139 : *
140 : * Now all rotations and mirroring can appear. This are resolved constructing
141 : * a mapping.
142 : */
143 :
144 0 : void MapLagrangeMultiIndexQuad(std::vector<size_t>& vOrientOffset,
145 : const int id0, bool sameOrientation,
146 : const size_t pOuter)
147 : {
148 : // in the inner, the number of dofs is as if it would be an element of order p-2.
149 0 : const size_t p = pOuter-2;
150 :
151 : // resize array
152 : vOrientOffset.clear();
153 0 : vOrientOffset.reserve((p+1)*(p+1));
154 :
155 : // loop mapped indices as required by rotated face
156 0 : for(size_t mj = 0; mj <= p; ++mj){
157 0 : for(size_t mi = 0; mi <= p; ++mi){
158 :
159 : // get corresponding multiindex in "natural" numbering
160 : size_t i,j;
161 0 : switch(id0)
162 : {
163 : case 0: i = mi; j = mj; break;
164 0 : case 1: i = mj; j = p-mi; break;
165 0 : case 2: i = p-mi; j = p-mj; break;
166 0 : case 3: i = p-mj; j = mi; break;
167 0 : default: UG_THROW("Orientation Quad: Corner "<<id0<<" invalid.");
168 : }
169 0 : if(!sameOrientation) std::swap(i, j);
170 :
171 : // linearize index
172 0 : const size_t naturalIndex = (p+1) * j + i;
173 :
174 : // set mapping
175 0 : vOrientOffset.push_back(naturalIndex);
176 : }
177 : }
178 0 : };
179 :
180 0 : void MapLagrangeMultiIndexTriangle(std::vector<size_t>& vOrientOffset,
181 : const int id0, bool sameOrientation,
182 : const size_t pOuter)
183 : {
184 : // in the inner, the number of dofs is as if it would be an element of order p-3.
185 0 : const size_t p = pOuter-3;
186 :
187 : // resize array
188 : vOrientOffset.clear();
189 :
190 : // loop mapped indices as required by rotated face
191 0 : for(size_t mj = 0; mj <= p; ++mj){
192 0 : for(size_t mi = 0; mi <= p-mj; ++mi){
193 :
194 : // get corresponding multiindex in "natural" numbering
195 : size_t i,j;
196 0 : switch(id0)
197 : {
198 : case 0: i = mi; j = mj; break;
199 0 : case 1: i = mj; j = p-mi-mj; break;
200 0 : case 2: i = p-mi-mj; j = mi; break;
201 0 : default: UG_THROW("Orientation Triangle: Corner "<<id0<<" invalid.");
202 : }
203 0 : if(!sameOrientation) std::swap(i, j);
204 :
205 : // linearize index and mapped index
206 0 : size_t naturalIndex = i;
207 0 : for(size_t c = 0; c < j; ++c)
208 0 : naturalIndex += (p+1-c);
209 :
210 : // set mapping
211 0 : vOrientOffset.push_back(naturalIndex);
212 : }
213 : }
214 0 : };
215 :
216 0 : void ComputeOrientationOffsetLagrange(std::vector<size_t>& vOrientOffset,
217 : Volume* volume, Face* face, size_t nrFace,
218 : const size_t p)
219 : {
220 : // get face descriptor
221 0 : FaceDescriptor fd;
222 0 : volume->face_desc(nrFace, fd);
223 :
224 : // find id0 and orientation
225 0 : const int numCo = face->num_vertices();
226 0 : const int id0 = GetVertexIndex(&fd, face->vertex(0));
227 0 : const int id1 = GetVertexIndex(&fd, face->vertex(1));
228 0 : const bool sameOrientation = (id1 == (id0+1)%numCo);
229 :
230 0 : switch(numCo){
231 0 : case 3:
232 0 : MapLagrangeMultiIndexTriangle(vOrientOffset, id0, sameOrientation, p);
233 : break;
234 0 : case 4:
235 0 : MapLagrangeMultiIndexQuad(vOrientOffset, id0, sameOrientation, p);
236 : break;
237 0 : default: UG_THROW("No corner number "<<numCo<<" implemented.");
238 : }
239 0 : };
240 :
241 0 : void ComputeOrientationOffsetLagrange(std::vector<size_t>& vOrientOffset,
242 : GridObject* volume, GridObject* face, size_t nrFace,
243 : const size_t p)
244 : {
245 0 : UG_THROW("Should never be called.")
246 : }
247 :
248 :
249 : ////////////////////////////////////////////////////////////////////////////////
250 : // General Implementation
251 : ////////////////////////////////////////////////////////////////////////////////
252 :
253 : template <typename TBaseElem, typename TSubBaseElem>
254 0 : void ComputeOrientationOffsetGeneric(std::vector<size_t>& vOrientOffset,
255 : TBaseElem* elem, TSubBaseElem* sub, size_t nrSub,
256 : const LFEID& lfeid)
257 : {
258 : vOrientOffset.clear();
259 :
260 : // if subelem higher dim than elem, no orientation
261 : if(TSubBaseElem::dim >= TBaseElem::dim) return;
262 :
263 0 : switch(lfeid.type()){
264 : // lagrange: orientate
265 : case LFEID::LAGRANGE:
266 : // only orientate if sub-dim < fct-space dim
267 0 : if(!(TSubBaseElem::dim < lfeid.dim())) return;
268 :
269 : // only orientate for p > 2, (else only max 1 DoF on sub)
270 0 : if(lfeid.order() <= 2) return;
271 :
272 : // orientate
273 0 : ComputeOrientationOffsetLagrange(vOrientOffset, elem, sub, nrSub, lfeid.order());
274 0 : break;
275 :
276 : // other cases: no orientation
277 : default: return;
278 : }
279 : }
280 :
281 :
282 0 : void ComputeOrientationOffset(std::vector<size_t>& vOrientOffset,
283 : Volume* vol, Face* face, size_t nrFace,
284 : const LFEID& lfeid)
285 : {
286 0 : ComputeOrientationOffsetGeneric(vOrientOffset, vol, face, nrFace, lfeid);
287 0 : }
288 :
289 0 : void ComputeOrientationOffset(std::vector<size_t>& vOrientOffset,
290 : Volume* vol, Edge* edge, size_t nrEdge,
291 : const LFEID& lfeid)
292 : {
293 0 : ComputeOrientationOffsetGeneric(vOrientOffset, vol, edge, nrEdge, lfeid);
294 0 : }
295 :
296 0 : void ComputeOrientationOffset(std::vector<size_t>& vOrientOffset,
297 : Face* face, Edge* edge, size_t nrEdge,
298 : const LFEID& lfeid)
299 : {
300 0 : ComputeOrientationOffsetGeneric(vOrientOffset, face, edge, nrEdge, lfeid);
301 0 : }
302 :
303 0 : void ComputeOrientationOffset(std::vector<size_t>& vOrientOffset,
304 : GridObject* Elem, GridObject* SubElem, size_t nrSub,
305 : const LFEID& lfeid)
306 : {
307 : // general case: no offset needed
308 : vOrientOffset.clear();
309 0 : }
310 :
311 : } // end namespace ug
|
