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 : #ifndef __H__UG__LIB_DISC__COMMON__GEOMETRY_UTIL__
34 : #define __H__UG__LIB_DISC__COMMON__GEOMETRY_UTIL__
35 :
36 : #include <vector>
37 : #include <cmath>
38 :
39 : #include "common/common.h"
40 : #include "lib_disc/reference_element/reference_element.h"
41 : #include "lib_disc/reference_element/element_list_traits.h"
42 : #include "lib_disc/domain_traits.h"
43 : #include "common/util/provider.h"
44 :
45 : namespace ug{
46 :
47 : ///////////////////////////////////////////////////////////////
48 : /// Volume of an Element in a given Dimension
49 : template <typename TRefElem, int TWorldDim>
50 : inline number ElementSize(const MathVector<TWorldDim>* vCornerCoords);
51 :
52 : ///////////////////////////////////////////////////////////////
53 : ///////////////////////////////////////////////////////////////
54 : // 1D Reference Elements
55 : ///////////////////////////////////////////////////////////////
56 : ///////////////////////////////////////////////////////////////
57 :
58 : ///////////////////////////////////////////////////////////////
59 : /// Volume of a Line in 1d
60 : /**
61 : * This function returns the volume of a line in 1d.
62 : *
63 : * \param[in] vCornerCoords Vector of corner coordinates (2 corners)
64 : * \return number Volume of Line
65 : */
66 : template <>
67 : inline number ElementSize<ReferenceEdge, 1>(const MathVector<1>* vCornerCoords)
68 : {
69 : return VecDistance(vCornerCoords[0], vCornerCoords[1]);
70 : }
71 :
72 : ///////////////////////////////////////////////////////////////
73 : /// Volume of a Line in 2d
74 : /**
75 : * This function returns the volume of a line in 2d.
76 : *
77 : * \param[in] vCornerCoords Vector of corner coordinates (2 corners)
78 : * \return number Volume of Line
79 : */
80 : template <>
81 : inline number ElementSize<ReferenceEdge, 2>(const MathVector<2>* vCornerCoords)
82 : {
83 0 : return VecDistance(vCornerCoords[0], vCornerCoords[1]);
84 : }
85 :
86 : ///////////////////////////////////////////////////////////////
87 : /// Volume of a Line in 3d
88 : /**
89 : * This function returns the volume of a line in 3d.
90 : *
91 : * \param[in] vCornerCoords Vector of corner coordinates (2 corners)
92 : * \return number Volume of Line
93 : */
94 : template <>
95 0 : inline number ElementSize<ReferenceEdge, 3>(const MathVector<3>* vCornerCoords)
96 : {
97 0 : return VecDistance(vCornerCoords[0], vCornerCoords[1]);
98 : }
99 :
100 : ///////////////////////////////////////////////////////////////
101 : ///////////////////////////////////////////////////////////////
102 : // 2D Reference Elements
103 : ///////////////////////////////////////////////////////////////
104 : ///////////////////////////////////////////////////////////////
105 :
106 : ///////////////////////////////////////////////////////////////
107 : /// Volume of a Triangle in 2d
108 : /**
109 : * This function returns the volume of a triangle in 2d.
110 : * The Volume is computed via:
111 : * F = 1/2 * fabs( (x2-x1)*(y1-y0) - (x1-x0)*(y2-y0) )
112 : *
113 : *
114 : * \param[in] vCornerCoords Vector of corner coordinates (3 corners)
115 : * \return number Volume of Triangle
116 : */
117 : template <>
118 : inline number ElementSize<ReferenceTriangle, 2>(const MathVector<2>* vCornerCoords)
119 : {
120 0 : return(0.5*fabs((vCornerCoords[1][1]-vCornerCoords[0][1])*(vCornerCoords[2][0]-vCornerCoords[0][0])
121 0 : -(vCornerCoords[1][0]-vCornerCoords[0][0])*(vCornerCoords[2][1]-vCornerCoords[0][1])));
122 : }
123 :
124 : ///////////////////////////////////////////////////////////////
125 : /// Volume of a Triangle in 3d
126 : /**
127 : * This function returns the volume of a triangle in 3d.
128 : * The Volume is computed via:
129 : * F = 1/2 * | (c - a) x (b - a) | (Here, a,b,c are the vectors of the corners)
130 : *
131 : * \param[in] vCornerCoords Vector of corner coordinates (3 corners)
132 : * \return number Volume of Triangle
133 : */
134 : template <>
135 0 : inline number ElementSize<ReferenceTriangle, 3>(const MathVector<3>* vCornerCoords)
136 : {
137 : MathVector<3> x10, x20, n;
138 :
139 : VecSubtract(x10, vCornerCoords[1], vCornerCoords[0]);
140 : VecSubtract(x20, vCornerCoords[2], vCornerCoords[0]);
141 0 : VecCross(n, x10, x20);
142 :
143 0 : return(0.5 * VecTwoNorm(n) );
144 : }
145 :
146 : ///////////////////////////////////////////////////////////////
147 : /// Volume of a Quadrilateral in 2d
148 : /**
149 : * This function returns the volume of a quadrilateral in 2d.
150 : * The Volume is computed via:
151 : * F = 1/2 * | (y3-y1)*(x2-x0) - (x3-x1)*(y2-y0) |
152 : *
153 : * \param[in] vCornerCoords Vector of corner coordinates (4 corners)
154 : * \return number Volume of Quadrilateral
155 : */
156 : template <>
157 : inline number ElementSize<ReferenceQuadrilateral, 2>(const MathVector<2>* vCornerCoords)
158 : {
159 0 : return( 0.5*fabs( (vCornerCoords[3][1]-vCornerCoords[1][1])*(vCornerCoords[2][0]-vCornerCoords[0][0])
160 0 : -(vCornerCoords[3][0]-vCornerCoords[1][0])*(vCornerCoords[2][1]-vCornerCoords[0][1]) ) );
161 : }
162 :
163 : ///////////////////////////////////////////////////////////////
164 : /// Volume of a Quadrilateral in 3d
165 : /**
166 : * This function returns the volume of a quadrilateral in 3d.
167 : * The Volume is computed via:
168 : * F = 1/2 * | (c - a) x (d - b) | (Here, a,b,c,d are the vectors of the corners)
169 : *
170 : * The corner coordinates must be given in counter - clockwise order
171 : *
172 : * \param[in] vCornerCoords Vector of corner coordinates (4 corners)
173 : * \return number Volume of Quadrilateral
174 : */
175 : template <>
176 0 : inline number ElementSize<ReferenceQuadrilateral, 3>(const MathVector<3>* vCornerCoords)
177 : {
178 : MathVector<3> x20, x31, n;
179 :
180 : VecSubtract(x20, vCornerCoords[2], vCornerCoords[0]);
181 : VecSubtract(x31, vCornerCoords[3], vCornerCoords[1]);
182 0 : VecCross(n, x20, x31);
183 :
184 0 : return(0.5 * VecTwoNorm(n) );
185 : }
186 :
187 :
188 : ///////////////////////////////////////////////////////////////
189 : ///////////////////////////////////////////////////////////////
190 : // 3D Reference Elements
191 : ///////////////////////////////////////////////////////////////
192 : ///////////////////////////////////////////////////////////////
193 :
194 : ///////////////////////////////////////////////////////////////
195 : /// Volume of a Tetrahedron in 3d
196 : /**
197 : * This function returns the volume of a tetrahedron in 3d.
198 : * The Volume is computed via:
199 : * V = 1/6 * | ( (b - a) x (c - a) ) * (d - a) | (Here, a,b,c,d are the vectors of the corners)
200 : *
201 : * This is motivated by the formula: V = 1/3 * (S * h) with
202 : * - S is the area of the base
203 : * - h is the height
204 : *
205 : * \param[in] vCornerCoords Vector of corner coordinates (4 corners)
206 : * \return number Volume of Tetrahedron
207 : */
208 : template <>
209 0 : inline number ElementSize<ReferenceTetrahedron, 3>(const MathVector<3>* vCornerCoords)
210 : {
211 : MathVector<3> x10, x20, x30, n;
212 :
213 : VecSubtract(x10, vCornerCoords[1], vCornerCoords[0]);
214 : VecSubtract(x20, vCornerCoords[2], vCornerCoords[0]);
215 : VecSubtract(x30, vCornerCoords[3], vCornerCoords[0]);
216 0 : VecCross(n, x10, x20);
217 :
218 0 : return (1./6.) * fabs ( VecDot(n, x30) );
219 : }
220 :
221 : ///////////////////////////////////////////////////////////////
222 : /// Volume of a Pyramid in 3d
223 : /**
224 : * This function returns the volume of a pyramid in 3d.
225 : * The volume is computed via: V = 1/3 * (S * h) with
226 : * - S is the area of the base
227 : * - h is the height
228 : *
229 : * The corner coordinates must be given as prescribed by the reference element
230 : *
231 : * \param[in] vCornerCoords Vector of corner coordinates (5 corners)
232 : * \return number Volume of Pyramid
233 : */
234 : template <>
235 0 : inline number ElementSize<ReferencePyramid, 3>(const MathVector<3>* vCornerCoords)
236 : {
237 : MathVector<3> x20, x31, x40, n;
238 :
239 : VecSubtract(x20, vCornerCoords[2], vCornerCoords[0]);
240 : VecSubtract(x31, vCornerCoords[3], vCornerCoords[1]);
241 : VecSubtract(x40, vCornerCoords[4], vCornerCoords[0]);
242 0 : VecCross(n, x20, x31);
243 :
244 0 : return (1./6.) * VecDot(n, x40);
245 : }
246 :
247 : ///////////////////////////////////////////////////////////////
248 : /// Volume of a Prism in 3d
249 : /**
250 : * This function returns the volume of a prism in 3d. Therefore, the
251 : * element is divided into pyramid {x0, x_1, x_4, x_3; x_5} and a
252 : * tetrahedron {x_0, x_1, x_2; x_5}, whose volumes are computed and added.
253 : *
254 : * The corner coordinates must be given as prescribed by the reference element
255 : *
256 : * \param[in] vCornerCoords Vector of corner coordinates (6 corners)
257 : * \return number Volume of Prism
258 : */
259 : template <>
260 0 : inline number ElementSize<ReferencePrism, 3>(const MathVector<3>* vCornerCoords)
261 : {
262 : MathVector<3> x40, x13, x10, x20, x50, m, n;
263 :
264 : VecSubtract(x40, vCornerCoords[4], vCornerCoords[0]);
265 : VecSubtract(x13, vCornerCoords[1], vCornerCoords[3]);
266 : VecSubtract(x10, vCornerCoords[1], vCornerCoords[0]);
267 : VecSubtract(x20, vCornerCoords[2], vCornerCoords[0]);
268 :
269 : // height for both subelements
270 : VecSubtract(x50, vCornerCoords[5], vCornerCoords[0]);
271 :
272 0 : VecCross(m, x40, x13); // base of pyramid
273 0 : VecCross(n, x10, x20); // base of tetrahedron
274 :
275 : // n = n + m
276 : VecAppend(n, m);
277 :
278 0 : return (1./6.) * VecDot(n, x50);
279 : }
280 :
281 : ///////////////////////////////////////////////////////////////
282 : /// Volume of a Hexahedron in 3d
283 : /**
284 : * This function returns the volume of a hexahedron in 3d. Therefore, the
285 : * element is divided into two prisms:
286 : * 1) {x0, x_1, x_2, x_4; x_5, x_6} and
287 : * 2) {x0, x_2, x_3, x_4; x_6, x_7}
288 : *
289 : * The corner coordinates must be given as prescribed by the reference element
290 : *
291 : * \param[in] vCornerCoords Vector of corner coordinates (8 corners)
292 : * \return number Volume of Hexahedron
293 : */
294 : template <>
295 0 : inline number ElementSize<ReferenceHexahedron, 3>(const MathVector<3>* vCornerCoords)
296 : {
297 : MathVector<3> x50, x14, x10, x20, x60, m1, n1;
298 : MathVector<3> x24, x30, x70, m2, n2;
299 :
300 : // 1. prism
301 : VecSubtract(x50, vCornerCoords[5], vCornerCoords[0]);
302 : VecSubtract(x14, vCornerCoords[1], vCornerCoords[4]);
303 : VecSubtract(x10, vCornerCoords[1], vCornerCoords[0]);
304 : VecSubtract(x20, vCornerCoords[2], vCornerCoords[0]);
305 : VecSubtract(x60, vCornerCoords[6], vCornerCoords[0]);
306 :
307 0 : VecCross(m1, x50, x14); // base of pyramid
308 0 : VecCross(n1, x10, x20); // base of tetrahedron
309 :
310 : // n1 = n1 + m1
311 : VecAppend(n1, m1);
312 :
313 : // 2. prism
314 : //VecSubtract(x60, vCornerCoords[6], vCornerCoords[0]);
315 : VecSubtract(x24, vCornerCoords[2], vCornerCoords[4]);
316 : //VecSubtract(x20, vCornerCoords[2], vCornerCoords[0]);
317 : VecSubtract(x30, vCornerCoords[3], vCornerCoords[0]);
318 : VecSubtract(x70, vCornerCoords[7], vCornerCoords[0]);
319 :
320 0 : VecCross(m2, x60, x24); // base of pyramid
321 0 : VecCross(n2, x20, x30); // base of tetrahedron
322 :
323 : // n2 = n2 + m2
324 : VecAppend(n2, m2);
325 :
326 0 : return (1./6.) * (VecDot(n1, x60) + VecDot(n2, x70));
327 : }
328 :
329 : ///////////////////////////////////////////////////////////////
330 : /// Volume of an Octahedron in 3d
331 : /**
332 : * This function returns the volume of an octhedron in 3d
333 : * by calculating the volumes of the upper and lower pyramid
334 : * the octahedron consists of.
335 : * The pyramidal volumes are computed via: V = 1/3 * (S * h) with
336 : * - S is the area of the base
337 : * - h is the height
338 : *
339 : * The corner coordinates must be given as prescribed by the reference element
340 : *
341 : * \param[in] vCornerCoords Vector of corner coordinates (6 corners)
342 : * \return number Volume of Octahedron
343 : */
344 : template <>
345 0 : inline number ElementSize<ReferenceOctahedron, 3>(const MathVector<3>* vCornerCoords)
346 : {
347 : MathVector<3> x31, x42, x51, n;
348 : MathVector<3> x01;
349 :
350 : VecSubtract(x31, vCornerCoords[3], vCornerCoords[1]);
351 : VecSubtract(x42, vCornerCoords[4], vCornerCoords[2]);
352 : VecSubtract(x51, vCornerCoords[5], vCornerCoords[1]);
353 0 : VecCross(n, x31, x42);
354 :
355 : //VecSubtract(x31, vCornerCoords[3], vCornerCoords[1]);
356 : //VecSubtract(x42, vCornerCoords[4], vCornerCoords[2]);
357 : VecSubtract(x01, vCornerCoords[0], vCornerCoords[1]);
358 : //VecCross(n, x31, x42);
359 :
360 0 : number volTopPyr = (1./6.) * fabs(VecDot(n, x51));
361 0 : number volBottomPyr = (1./6.) * fabs(VecDot(n, x01));
362 :
363 0 : return volTopPyr + volBottomPyr;
364 : }
365 :
366 : ///////////////////////////////////////////////////////////////
367 : // run-time size of element
368 : ///////////////////////////////////////////////////////////////
369 :
370 : template <int dim>
371 : inline number ElementSize(ReferenceObjectID roid, const MathVector<dim>* vCornerCoords);
372 :
373 : template <>
374 0 : inline number ElementSize<1>(ReferenceObjectID roid, const MathVector<1>* vCornerCoords)
375 : {
376 0 : switch(roid)
377 : {
378 : case ROID_VERTEX: return 1.0;
379 0 : case ROID_EDGE: return ElementSize<ReferenceEdge, 1>(vCornerCoords);
380 0 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 1.");
381 : }
382 : }
383 :
384 : template <>
385 0 : inline number ElementSize<2>(ReferenceObjectID roid, const MathVector<2>* vCornerCoords)
386 : {
387 0 : switch(roid)
388 : {
389 : case ROID_VERTEX: return 1.0;
390 0 : case ROID_EDGE: return ElementSize<ReferenceEdge, 2>(vCornerCoords);
391 0 : case ROID_TRIANGLE: return ElementSize<ReferenceTriangle, 2>(vCornerCoords);
392 0 : case ROID_QUADRILATERAL: return ElementSize<ReferenceQuadrilateral, 2>(vCornerCoords);
393 0 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 2.");
394 : }
395 : }
396 :
397 : template <>
398 0 : inline number ElementSize<3>(ReferenceObjectID roid, const MathVector<3>* vCornerCoords)
399 : {
400 0 : switch(roid)
401 : {
402 : case ROID_VERTEX: return 1.0;
403 0 : case ROID_EDGE: return ElementSize<ReferenceEdge, 3>(vCornerCoords);
404 0 : case ROID_TRIANGLE: return ElementSize<ReferenceTriangle, 3>(vCornerCoords);
405 0 : case ROID_QUADRILATERAL: return ElementSize<ReferenceQuadrilateral, 3>(vCornerCoords);
406 0 : case ROID_TETRAHEDRON: return ElementSize<ReferenceTetrahedron, 3>(vCornerCoords);
407 0 : case ROID_PYRAMID: return ElementSize<ReferencePyramid, 3>(vCornerCoords);
408 0 : case ROID_PRISM: return ElementSize<ReferencePrism, 3>(vCornerCoords);
409 0 : case ROID_HEXAHEDRON: return ElementSize<ReferenceHexahedron, 3>(vCornerCoords);
410 0 : case ROID_OCTAHEDRON: return ElementSize<ReferenceOctahedron, 3>(vCornerCoords);
411 0 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 3.");
412 : }
413 : }
414 :
415 : ///////////////////////////////////////////////////////////////
416 : ///////////////////////////////////////////////////////////////
417 : // Normals on Elements
418 : ///////////////////////////////////////////////////////////////
419 : ///////////////////////////////////////////////////////////////
420 :
421 : ///////////////////////////////////////////////////////////////
422 : /// Normal to an Element in a given Dimension
423 : template <typename TRefElem, int TWorldDim>
424 : inline void ElementNormal(MathVector<TWorldDim>& normalOut, const MathVector<TWorldDim>* vCornerCoords);
425 :
426 : ///////////////////////////////////////////////////////////////
427 : /// Normal to a Point in 1d
428 : /**
429 : * This function returns the normal of a point in 1d.
430 : * This always 1.0.
431 : *
432 : * \param[in] vCornerCoords Vector of corner coordinates
433 : * \param[out] normalOut Normal
434 : */
435 : template <>
436 : inline void ElementNormal<ReferenceVertex, 1>(MathVector<1>& normalOut, const MathVector<1>* vCornerCoords)
437 : {
438 0 : normalOut[0] = 1.0;
439 : }
440 :
441 : ///////////////////////////////////////////////////////////////
442 : /// Normal to a Point in 2d
443 : /**
444 : * This function returns the normal of a point in 2d.
445 : * This can only be understood as the normal on a corner of an edge in 2d.
446 : * Therefore, it will be assumed that vCornerCoords has exactly 2 entries
447 : * defining the edge. Then the normal is the one-dim. normal on the first vertex
448 : * embedded in the 1D subspace defined by the edge.
449 : *
450 : * \param[in] vCornerCoords Vector of corner coordinates
451 : * \param[out] normalOut Normal
452 : */
453 : template <>
454 0 : inline void ElementNormal<ReferenceVertex, 2>(MathVector<2>& normalOut, const MathVector<2>* vCornerCoords)
455 : {
456 : try
457 : {
458 : normalOut = vCornerCoords[0];
459 : normalOut -= vCornerCoords[1];
460 0 : VecNormalize(normalOut, normalOut);
461 : }
462 : UG_CATCH_THROW("Element normal can not be computed. Maybe there are not enough vertices to work on.")
463 0 : }
464 :
465 : ///////////////////////////////////////////////////////////////
466 : /// Normal to a Point in 3d
467 : /**
468 : * This function returns the normal of a point in 3d.
469 : * This can only be understood as the normal on a corner of an edge in 3d.
470 : * Therefore, it will be assumed that vCornerCoords has exactly 2 entries
471 : * defining the edge. Then the normal is the one-dim. normal on the first vertex
472 : * embedded in the 1D subspace defined by the edge.
473 : *
474 : * \param[in] vCornerCoords Vector of corner coordinates
475 : * \param[out] normalOut Normal
476 : */
477 : template <>
478 0 : inline void ElementNormal<ReferenceVertex, 3>(MathVector<3>& normalOut, const MathVector<3>* vCornerCoords)
479 : {
480 : try
481 : {
482 : normalOut = vCornerCoords[0];
483 : normalOut -= vCornerCoords[1];
484 0 : VecNormalize(normalOut, normalOut);
485 : }
486 : UG_CATCH_THROW("Element normal can not be computed. Maybe there are not enough vertices to work on.")
487 0 : }
488 :
489 : ///////////////////////////////////////////////////////////////
490 : /// Normal to a Line in 2d
491 : /**
492 : * This function returns the normal of a line in 2d.
493 : * The normal is in clockwise direction to the vector (x1-x0).
494 : * The norm of the normal is the size of the line.
495 : *
496 : * \param[in] vCornerCoords Vector of corner coordinates
497 : * \param[out] normalOut Normal
498 : */
499 : template <>
500 : inline void ElementNormal<ReferenceEdge, 2>(MathVector<2>& normalOut, const MathVector<2>* vCornerCoords)
501 : {
502 : MathVector<2> diff(vCornerCoords[1]);
503 : diff -= vCornerCoords[0];
504 :
505 0 : normalOut[0] = diff[1];
506 0 : normalOut[1] = -diff[0];
507 0 : }
508 :
509 : ///////////////////////////////////////////////////////////////
510 : /// Normal to a Line in 3d
511 : /**
512 : * This function returns the normal of a line in 3d.
513 : * This can only be understood as the normal on an edge in a 2d manifold
514 : * defined by at least a triangle. Therefore, it is assumed that vCornerCoords
515 : * contains at least three vertices.
516 : * The normal is computed as the outer normal on the first edge (v2-v1) of the triangle
517 : * and embedded in the 2d subspace defined by the triangle.
518 : * The norm of the normal is the size of the line.
519 : *
520 : * \param[in] vCornerCoords Vector of corner coordinates
521 : * \param[out] normalOut Normal
522 : */
523 : template <>
524 0 : inline void ElementNormal<ReferenceEdge, 3>(MathVector<3>& normalOut, const MathVector<3>* vCornerCoords)
525 : {
526 : try
527 : {
528 : // normal on triangle
529 : MathVector<3> edge0, edge1;
530 : VecSubtract(edge0, vCornerCoords[1], vCornerCoords[0]);
531 : VecSubtract(edge1, vCornerCoords[2], vCornerCoords[1]);
532 0 : VecCross(normalOut, edge0, edge1);
533 : // normal an edge is edge x normal on triangle
534 0 : VecCross(normalOut, edge0, normalOut);
535 : // scale
536 0 : VecNormalize(normalOut, normalOut);
537 : VecScale(normalOut, normalOut, VecTwoNorm(edge0));
538 : }
539 : UG_CATCH_THROW("Element normal can not be computed. Maybe there are not enough vertices to work on.")
540 0 : }
541 :
542 : ///////////////////////////////////////////////////////////////
543 : /// Normal to a Triangle in 3d
544 : /**
545 : * This function returns the normal of a triangle in 3d.
546 : * The orientation is right-handed (i.e. forming with the four fingers of the
547 : * right hand a circle in direction of the corner nodes, the thumb points in
548 : * the direction of the normal)
549 : *
550 : * The norm of the normal is the size of the triangle.
551 : *
552 : * \param[in] vCornerCoords Vector of corner coordinates
553 : * \param[out] normalOut Normal
554 : */
555 : template <>
556 0 : inline void ElementNormal<ReferenceTriangle, 3>(MathVector<3>& normalOut, const MathVector<3>* vCornerCoords)
557 : {
558 : MathVector<3> a, b;
559 : VecSubtract(a, vCornerCoords[1], vCornerCoords[0]);
560 : VecSubtract(b, vCornerCoords[2], vCornerCoords[0]);
561 0 : VecCross(normalOut, a,b);
562 : VecScale(normalOut, normalOut, 0.5);
563 0 : }
564 :
565 : ///////////////////////////////////////////////////////////////
566 : /// Normal to a Quadrilateral in 3d
567 : /**
568 : * This function returns the normal of a quadrilateral in 3d.
569 : * The orientation is right-handed (i.e. forming with the four fingers of the
570 : * right hand a circle in direction of the corner nodes, the thumb points in
571 : * the direction of the normal)
572 : *
573 : * The norm of the normal is the size of the quadrilateral.
574 : *
575 : * \param[in] vCornerCoords Vector of corner coordinates
576 : * \param[out] normalOut Normal
577 : */
578 : template <>
579 0 : inline void ElementNormal<ReferenceQuadrilateral, 3>(MathVector<3>& normalOut, const MathVector<3>* vCornerCoords)
580 : {
581 : MathVector<3> a, b;
582 : VecSubtract(a, vCornerCoords[2], vCornerCoords[0]);
583 : VecSubtract(b, vCornerCoords[3], vCornerCoords[1]);
584 0 : VecCross(normalOut, a,b);
585 : VecScale(normalOut, normalOut, 0.5);
586 0 : }
587 :
588 : ///////////////////////////////////////////////////////////////
589 : // run-time normal of element
590 : ///////////////////////////////////////////////////////////////
591 :
592 : template <int dim>
593 : inline void ElementNormal(ReferenceObjectID roid, MathVector<dim>& normalOut, const MathVector<dim>* vCornerCoords);
594 :
595 : template <>
596 0 : inline void ElementNormal<1>(ReferenceObjectID roid, MathVector<1>& normalOut, const MathVector<1>* vCornerCoords)
597 : {
598 0 : switch(roid)
599 : {
600 0 : case ROID_VERTEX: ElementNormal<ReferenceVertex, 1>(normalOut, vCornerCoords); return;
601 0 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 1.");
602 : }
603 : }
604 :
605 : template <>
606 0 : inline void ElementNormal<2>(ReferenceObjectID roid, MathVector<2>& normalOut, const MathVector<2>* vCornerCoords)
607 : {
608 0 : switch(roid)
609 : {
610 0 : case ROID_VERTEX: ElementNormal<ReferenceVertex, 2>(normalOut, vCornerCoords); return;
611 : case ROID_EDGE: ElementNormal<ReferenceEdge, 2>(normalOut, vCornerCoords); return;
612 0 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 2.");
613 : }
614 : }
615 :
616 : template <>
617 0 : inline void ElementNormal<3>(ReferenceObjectID roid, MathVector<3>& normalOut, const MathVector<3>* vCornerCoords)
618 : {
619 0 : switch(roid)
620 : {
621 0 : case ROID_VERTEX: ElementNormal<ReferenceVertex, 3>(normalOut, vCornerCoords); return;
622 0 : case ROID_EDGE: ElementNormal<ReferenceEdge, 3>(normalOut, vCornerCoords); return;
623 0 : case ROID_TRIANGLE: ElementNormal<ReferenceTriangle, 3>(normalOut, vCornerCoords); return;
624 0 : case ROID_QUADRILATERAL: ElementNormal<ReferenceQuadrilateral, 3>(normalOut, vCornerCoords); return;
625 0 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 3.");
626 : }
627 : }
628 :
629 : ///////////////////////////////////////////////////////////////
630 : /// Normal to a side of an Element in a given Dimension
631 : /**
632 : * This function computes the outer normal to a given side
633 : * of an element. The Euclidean norm of the normal is the
634 : * area of the side. Note that the normal is computed in
635 : * the same dimensionality as the reference element itself.
636 : *
637 : * \param[in] side index of the side of the element
638 : * \param[in] vCornerCoords array of the global coordinates of the corners
639 : * \param[out] normalOut the computed normal
640 : */
641 : template <typename TRefElem, int TWorldDim>
642 0 : inline void SideNormal(MathVector<TWorldDim>& normalOut, int side, const MathVector<TWorldDim>* vCornerCoords)
643 : {
644 : static const int dim = (int) TRefElem::dim;
645 : static const int maxSideCorners = element_list_traits<typename domain_traits<dim-1>::DimElemList>::maxCorners;
646 :
647 : // Get own reference element and the side roid:
648 : TRefElem & rRefElem = (TRefElem&) ReferenceElementProvider::get(TRefElem::REFERENCE_OBJECT_ID);
649 0 : ReferenceObjectID sideRoid = rRefElem.roid(dim-1,side);
650 :
651 : // Get the coordinates of the vertices:
652 0 : MathVector<TWorldDim> vSideCorner [(dim == TWorldDim)? maxSideCorners : maxSideCorners + 1];
653 : size_t numSideCorners = rRefElem.num(dim-1, side, 0);
654 0 : for (size_t co = 0; co < numSideCorners; ++co)
655 0 : vSideCorner[co] = vCornerCoords[rRefElem.id(dim-1, side, 0, co)];
656 : // we need another point if dim != TWorldDim:
657 : // take the highest-numbered corner of the next side, since
658 : // it is always different from the other points (is it not?)
659 : if (dim != TWorldDim)
660 : {
661 : vSideCorner[numSideCorners] =
662 0 : vCornerCoords[rRefElem.id(dim-1, (side+1)%rRefElem.num(dim-1), 0,
663 0 : rRefElem.num(dim-1, (side+1)%rRefElem.num(dim-1), 0)-1)];
664 : }
665 :
666 : // Get the normal:
667 0 : ElementNormal<TWorldDim>(sideRoid, normalOut, vSideCorner);
668 : // Note: We assume that the for the standard ordering, the last line computes
669 : // the outer normal.
670 0 : }
671 :
672 : /// Computation of the side normal for a generic reference element:
673 : template <int dim>
674 : inline void SideNormal(ReferenceObjectID roid, MathVector<dim>& normalOut, int side, const MathVector<dim>* vCornerCoords);
675 :
676 : template <>
677 : inline void SideNormal<1>(ReferenceObjectID roid, MathVector<1>& normalOut, int side, const MathVector<1>* vCornerCoords)
678 : {
679 : switch(roid)
680 : {
681 : case ROID_EDGE: SideNormal<ReferenceEdge,1>(normalOut, side, vCornerCoords); return;
682 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 1.");
683 : }
684 : }
685 :
686 : template <>
687 : inline void SideNormal<2>(ReferenceObjectID roid, MathVector<2>& normalOut, int side, const MathVector<2>* vCornerCoords)
688 : {
689 : switch(roid)
690 : {
691 : case ROID_EDGE: SideNormal<ReferenceEdge,2>(normalOut, side, vCornerCoords); return;
692 : case ROID_TRIANGLE: SideNormal<ReferenceTriangle,2>(normalOut, side, vCornerCoords); return;
693 : case ROID_QUADRILATERAL: SideNormal<ReferenceQuadrilateral,2>(normalOut, side, vCornerCoords); return;
694 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 2.");
695 : }
696 : }
697 :
698 : template <>
699 : inline void SideNormal<3>(ReferenceObjectID roid, MathVector<3>& normalOut, int side, const MathVector<3>* vCornerCoords)
700 : {
701 : switch(roid)
702 : {
703 : case ROID_EDGE: SideNormal<ReferenceEdge,3>(normalOut, side, vCornerCoords); return;
704 : case ROID_TRIANGLE: SideNormal<ReferenceTriangle,3>(normalOut, side, vCornerCoords); return;
705 : case ROID_QUADRILATERAL: SideNormal<ReferenceQuadrilateral,3>(normalOut, side, vCornerCoords); return;
706 : case ROID_TETRAHEDRON: SideNormal<ReferenceTetrahedron,3>(normalOut, side, vCornerCoords); return;
707 : case ROID_PYRAMID: SideNormal<ReferencePyramid,3>(normalOut, side, vCornerCoords); return;
708 : case ROID_PRISM: SideNormal<ReferencePrism,3>(normalOut, side, vCornerCoords); return;
709 : case ROID_HEXAHEDRON: SideNormal<ReferenceHexahedron,3>(normalOut, side, vCornerCoords); return;
710 : case ROID_OCTAHEDRON: SideNormal<ReferenceOctahedron,3>(normalOut, side, vCornerCoords); return;
711 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 3.");
712 : }
713 : }
714 :
715 : ///////////////////////////////////////////////////////////////
716 : ///////////////////////////////////////////////////////////////
717 : // ElementSideRayIntersection
718 : ///////////////////////////////////////////////////////////////
719 : ///////////////////////////////////////////////////////////////
720 :
721 : // wrapper class to distinguish reference dimesion
722 : template <typename TRefElem, int TWorldDim, int TRefDim = TRefElem::dim>
723 : struct ElementSideRayIntersectionWrapper
724 : {
725 0 : static bool apply( size_t& sideOut,
726 : MathVector<TWorldDim>& GlobalIntersectionPointOut,
727 : MathVector<TRefElem::dim>& LocalIntersectionPoint,
728 : const MathVector<TWorldDim>& From, const MathVector<TWorldDim>& Direction,
729 : bool bPositiv, const MathVector<TWorldDim>* vCornerCoords)
730 : {
731 0 : UG_THROW("Not implemented.");
732 : return false;
733 : }
734 : };
735 :
736 : // specialization for 2d
737 : template <typename TRefElem>
738 : struct ElementSideRayIntersectionWrapper<TRefElem, 2, 2>
739 : {
740 0 : static bool apply( size_t& sideOut,
741 : MathVector<2>& GlobalIntersectionPointOut,
742 : MathVector<TRefElem::dim>& LocalIntersectionPoint,
743 : const MathVector<2>& From, const MathVector<2>& Direction,
744 : bool bPositiv, const MathVector<2>* vCornerCoords)
745 : {
746 0 : static const TRefElem& rRefElem = Provider<TRefElem>::get();
747 :
748 : // reference dimension
749 : const int dim = TRefElem::dim;
750 :
751 : // parameters
752 0 : number bc = 0., t = 0.;
753 : size_t p0 = 0, p1 = 0;
754 :
755 : // find side
756 0 : for(sideOut = 0; sideOut < rRefElem.num(dim-1); ++sideOut)
757 : {
758 : // get corners
759 0 : p0 = rRefElem.id(dim-1, sideOut, 0, 0);
760 0 : p1 = rRefElem.id(dim-1, sideOut, 0, 1);
761 :
762 : // if match: break
763 0 : if(RayLineIntersection2d( GlobalIntersectionPointOut, bc, t,
764 0 : vCornerCoords[p0], vCornerCoords[p1],
765 : From, Direction))
766 : {
767 0 : if(bPositiv && t >= 0.0) break;
768 0 : else if(!bPositiv && t <= 0.0) break;
769 : }
770 : }
771 : // if not found
772 0 : if(sideOut >= rRefElem.num(dim-1))
773 0 : UG_THROW("ElementSideRayIntersection: no cut side found.");
774 :
775 : // Compute local intersection
776 0 : VecScaleAdd(LocalIntersectionPoint, bc, rRefElem.corner(p1), 1.-bc, rRefElem.corner(p0));
777 :
778 : // true if found
779 0 : return true;
780 : }
781 : };
782 :
783 : // specialization for 3d
784 : template <typename TRefElem>
785 : struct ElementSideRayIntersectionWrapper<TRefElem, 3, 3>
786 : {
787 0 : static bool apply( size_t& sideOut,
788 : MathVector<3>& GlobalIntersectionPointOut,
789 : MathVector<TRefElem::dim>& LocalIntersectionPoint,
790 : const MathVector<3>& From, const MathVector<3>& Direction,
791 : bool bPositiv, const MathVector<3>* vCornerCoords)
792 : {
793 0 : static const TRefElem& rRefElem = Provider<TRefElem>::get();
794 :
795 : // reference dimension
796 : const int dim = TRefElem::dim;
797 :
798 : // parameters
799 0 : number bc0 = 0., bc1 = 0., t = 0.;
800 : size_t p0 = 0, p1 = 0, p2 = 0;
801 :
802 : // find side
803 0 : for(sideOut = 0; sideOut < rRefElem.num(dim-1); ++sideOut)
804 : {
805 : // get corners
806 0 : p0 = rRefElem.id(dim-1, sideOut, 0, 0);
807 0 : p1 = rRefElem.id(dim-1, sideOut, 0, 1);
808 0 : p2 = rRefElem.id(dim-1, sideOut, 0, 2);
809 :
810 : // if match: break
811 0 : if(RayTriangleIntersection( GlobalIntersectionPointOut, bc0, bc1, t,
812 0 : vCornerCoords[p0], vCornerCoords[p1], vCornerCoords[p2],
813 : From, Direction))
814 : {
815 0 : if(bPositiv && t >= 0.0) break;
816 0 : else if(!bPositiv && t <= 0.0) break;
817 : }
818 :
819 : // second triangle (only if 4 corners)
820 0 : if(rRefElem.num(dim-1, sideOut, 0) == 3) continue;
821 :
822 : // get corner number 4
823 0 : p1 = rRefElem.id(dim-1, sideOut, 0, 3);
824 :
825 : // if match: break
826 0 : if(RayTriangleIntersection( GlobalIntersectionPointOut, bc0, bc1, t,
827 0 : vCornerCoords[p0], vCornerCoords[p1], vCornerCoords[p2],
828 : From, Direction))
829 : {
830 0 : if(bPositiv && t >= 0.0) break;
831 0 : else if(!bPositiv && t <= 0.0) break;
832 : }
833 : }
834 :
835 : // if not found
836 0 : if(sideOut >= rRefElem.num(dim-1))
837 0 : UG_THROW("ElementSideRayIntersection: no cut side found.");
838 :
839 : // Compute local intersection
840 0 : VecScaleAdd(LocalIntersectionPoint,
841 0 : (1.-bc0-bc1), rRefElem.corner(p0),
842 : bc0, rRefElem.corner(p1),
843 : bc1, rRefElem.corner(p2));
844 :
845 : // true if found
846 0 : return true;
847 : }
848 : };
849 :
850 :
851 : ///////////////////////////////////////////////////////////////
852 : /// ElementSideRayIntersection
853 : /**
854 : * This function computes the side of element, that is intersected
855 : * by a Ray given in Parameter form by 'From + c * Direction'. The
856 : * intersection is choose to be at positive parameter if bPositiv == true,
857 : * else the intersection at negative parameter is choosen.
858 : * Global and local coordinates of the intersection points are returned
859 : * as well as the reference element number of the side.
860 : *
861 : * \param[in] From Point of Ray
862 : * \param[in] Direction Direction of Ray
863 : * \param[in] bPositiv Flag, whether to search in positiv of negative direction
864 : * \param[in] vCornerCoords Vector of corner coordinates
865 : * \param[out] sideOut side of intersection
866 : * \param[out] GlobalIntersectionPointOut Intersection Point (global)
867 : * \param[out] LocalIntersectionPoint Intersection Point (local)
868 : */
869 : template <typename TRefElem, int TWorldDim>
870 : bool ElementSideRayIntersection( size_t& sideOut,
871 : MathVector<TWorldDim>& GlobalIntersectionPointOut,
872 : MathVector<TRefElem::dim>& LocalIntersectionPoint,
873 : const MathVector<TWorldDim>& From, const MathVector<TWorldDim>& Direction,
874 : bool bPositiv, const MathVector<TWorldDim>* vCornerCoords)
875 : {
876 : UG_ASSERT(VecTwoNorm(Direction) > 0, "Direction must be non-zero vector.");
877 : return ElementSideRayIntersectionWrapper<TRefElem, TWorldDim>::
878 0 : apply(sideOut, GlobalIntersectionPointOut, LocalIntersectionPoint,
879 0 : From, Direction, bPositiv, vCornerCoords);
880 : }
881 :
882 :
883 : ///////////////////////////////////////////////////////////////
884 : ///////////////////////////////////////////////////////////////
885 : // ElementSideRayIntersection
886 : ///////////////////////////////////////////////////////////////
887 : ///////////////////////////////////////////////////////////////
888 :
889 : // wrapper class to distinguish reference dimesion
890 : template <int TDim, int TWorldDim>
891 : struct SCVFofSCVRayIntersectionWrapper
892 : {
893 : static bool apply( size_t& sideOut, number& bc,
894 : MathVector<TWorldDim>& GlobalIntersectionPointOut,
895 : MathVector<TDim>& LocalIntersectionPoint,
896 : const MathVector<TWorldDim>& From, const MathVector<TWorldDim>& Direction,
897 : bool bPositiv, const MathVector<TWorldDim>* vCornerCoords)
898 : {
899 : UG_THROW("Not implemented.");
900 : return false;
901 : }
902 : };
903 :
904 : // specialization for 2d
905 : template <>
906 : struct SCVFofSCVRayIntersectionWrapper<2, 2>
907 : {
908 : static bool apply( size_t& sideOut, number& bc,
909 : MathVector<2>& GlobalIntersectionPointOut,
910 : MathVector<2>& LocalIntersectionPoint,
911 : const MathVector<2>& From, const MathVector<2>& Direction,
912 : bool bPositiv, const MathVector<2>* vCornerCoords)
913 : {
914 : static const ReferenceQuadrilateral& rRefElem = Provider<ReferenceQuadrilateral>::get();
915 :
916 : // reference dimension
917 : static const int dim = 2;
918 :
919 : // parameters
920 : bc = 0.;
921 : number t = 0.;
922 : size_t p0 = 0, p1 = 0;
923 :
924 : // find side
925 : for(sideOut = 0; sideOut < rRefElem.num(0); ++sideOut)
926 : {
927 : // get corners
928 : p0 = rRefElem.id(dim-1, sideOut, 0, 0);
929 : p1 = rRefElem.id(dim-1, sideOut, 0, 1);
930 :
931 : // if match: break
932 : if(RayLineIntersection2d( GlobalIntersectionPointOut, bc, t,
933 : vCornerCoords[p0], vCornerCoords[p1],
934 : From, Direction))
935 : {
936 : // skip if one root-point scvf
937 : if(fabs(t) <= std::numeric_limits<number>::epsilon() * 10)
938 : continue;
939 :
940 : // upwind / downwind switch
941 : if(bPositiv && t >= 0.0) break;
942 : else if(!bPositiv && t <= 0.0) break;
943 : }
944 : }
945 :
946 : // if not found
947 : if(sideOut >= rRefElem.num(0))
948 : UG_THROW("Side not found.");
949 :
950 : // Compute local intersection
951 : VecScaleAdd(LocalIntersectionPoint, bc, rRefElem.corner(p1), 1.-bc, rRefElem.corner(p0));
952 :
953 : // parameter of intersection should loop always from center to edgeMidpoint
954 : // thus, for side 2 this is correct, but for side 1 we have to invert direction
955 : if(sideOut == 1) bc = 1. - bc;
956 :
957 : // true if found on scvf
958 : if(sideOut == 1 || sideOut == 2) return true;
959 : // false if found on element corner
960 : else return false;
961 : }
962 : };
963 :
964 : ///////////////////////////////////////////////////////////////
965 : // SCVFofSCVRayIntersection
966 : /**
967 : * This function computes if another scvf of a scv is intersected
968 : * by a Ray given in Parameter form by 'Root + c * Direction'. The
969 : * intersection is choose to be at positive parameter if bPositiv == true,
970 : * else the intersection at negative parameter is choosen.
971 : * Global and local coordinates of the intersection points are returned
972 : * as well as the local number of the side of the scv that matches the intersection.
973 : *
974 : * \param[in] Root Point of Ray
975 : * \param[in] Direction Direction of Ray
976 : * \param[in] bPositiv Flag, whether to search in positiv of negative direction
977 : * \param[in] vCornerCoords Vector of corner coordinates
978 : * \param[out] sideOut side of intersection
979 : * \param[out] bc line parameter in [0,1] indicating position
980 : * of intersection point on scvf in direction
981 : * center to edge
982 : * \param[out] GlobalIntersectionPointOut Intersection Point (global)
983 : * \param[out] LocalIntersectionPoint Intersection Point (local)
984 : * \returns true if intersected with another scvf of the scv, false else
985 : */
986 : template <int TDim, int TWorldDim>
987 : bool SCVFofSCVRayIntersection( size_t& sideOut, number& bc,
988 : MathVector<TWorldDim>& GlobalIntersectionPointOut,
989 : MathVector<TDim>& LocalIntersectionPoint,
990 : const MathVector<TWorldDim>& Root, const MathVector<TWorldDim>& Direction,
991 : bool bPositiv, const MathVector<TWorldDim>* vCornerCoords)
992 : {
993 : UG_ASSERT(VecTwoNorm(Direction) > 0, "Direction must be non-zero vector.");
994 : return SCVFofSCVRayIntersectionWrapper<TDim, TWorldDim>::
995 : apply(sideOut, bc, GlobalIntersectionPointOut, LocalIntersectionPoint,
996 : Root, Direction, bPositiv, vCornerCoords);
997 : }
998 :
999 :
1000 :
1001 : ///////////////////////////////////////////////////////////////
1002 : /// Extension of an element (in a given dimension)
1003 : ///////////////////////////////////////////////////////////////
1004 :
1005 :
1006 : template <typename TVector>
1007 : inline void ComputeElementExtensionsSqForEdges(const TVector* vCornerCoords, TVector &ext)
1008 : {
1009 : VecElemProd(ext, vCornerCoords[0], vCornerCoords[1]);
1010 : }
1011 :
1012 :
1013 : template <int TWorldDim, int ncorners>
1014 : inline void ComputeElementExtensionsSq(const MathVector<TWorldDim>* vCornerCoords, MathVector<TWorldDim> &ext)
1015 : {
1016 : // compute center
1017 : MathVector<TWorldDim> mid = 0.0;
1018 : for (int i=ncorners-1; i>=0; --i){
1019 : VecAppend(mid, vCornerCoords[i]);
1020 : }
1021 : VecScale(mid, mid, 1.0/ncorners);
1022 :
1023 : // compute
1024 : ext = 0.0;
1025 : MathVector<TWorldDim> aux;
1026 : for (int i=ncorners-1; i>=0; --i){
1027 : VecSubtract(aux, mid, vCornerCoords[i]);
1028 : VecElemProd(aux, aux, aux);
1029 : VecAppend(ext, aux);
1030 : }
1031 : VecScale(ext, ext, 1.0/ncorners);
1032 : }
1033 :
1034 : ///////////////////////////////////////////////////////////////
1035 : // extensions of an element
1036 : ///////////////////////////////////////////////////////////////
1037 :
1038 :
1039 : template <int dim>
1040 : inline void ElementExtensionsSq(ReferenceObjectID roid, MathVector<dim>& ext, const MathVector<dim>* vCornerCoords);
1041 :
1042 : template <>
1043 : inline void ElementExtensionsSq<1>(ReferenceObjectID roid, MathVector<1>& ext, const MathVector<1>* vCornerCoords)
1044 : {
1045 : switch(roid)
1046 : {
1047 : case ROID_VERTEX: ext = 0; return;
1048 : case ROID_EDGE: ComputeElementExtensionsSqForEdges(vCornerCoords,ext); return;
1049 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 1.");
1050 : }
1051 : }
1052 :
1053 : template <>
1054 : inline void ElementExtensionsSq<2>(ReferenceObjectID roid, MathVector<2>& ext, const MathVector<2>* vCornerCoords)
1055 : {
1056 : switch(roid)
1057 : {
1058 : case ROID_VERTEX: ext = 0; return;
1059 : case ROID_EDGE: ComputeElementExtensionsSqForEdges(vCornerCoords,ext); return;
1060 : case ROID_TRIANGLE: ComputeElementExtensionsSq<2,3>(vCornerCoords,ext); return;
1061 : case ROID_QUADRILATERAL: ComputeElementExtensionsSq<2,4>(vCornerCoords,ext); return;
1062 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 2.");
1063 : }
1064 : }
1065 :
1066 : template <>
1067 : inline void ElementExtensionsSq<3>(ReferenceObjectID roid, MathVector<3>& ext, const MathVector<3>* vCornerCoords)
1068 : {
1069 : switch(roid)
1070 : {
1071 : case ROID_VERTEX: ext = 0; return;
1072 : case ROID_EDGE: ComputeElementExtensionsSqForEdges(vCornerCoords,ext); return;
1073 : case ROID_TRIANGLE: ComputeElementExtensionsSq<3,3>(vCornerCoords,ext); return;
1074 : case ROID_QUADRILATERAL: ComputeElementExtensionsSq<3,4>(vCornerCoords,ext); return;
1075 : case ROID_TETRAHEDRON: ComputeElementExtensionsSq<3,4>(vCornerCoords,ext); return;
1076 : case ROID_PYRAMID: ComputeElementExtensionsSq<3,5>(vCornerCoords,ext); return;
1077 : case ROID_PRISM: ComputeElementExtensionsSq<3,6>(vCornerCoords,ext); return;
1078 : case ROID_HEXAHEDRON: ComputeElementExtensionsSq<3,8>(vCornerCoords,ext); return;
1079 : case ROID_OCTAHEDRON: ComputeElementExtensionsSq<3,6>(vCornerCoords,ext); return;
1080 : default: UG_THROW("ReferenceObject "<<roid<<" not found in dim 3.");
1081 : }
1082 : }
1083 :
1084 :
1085 : } // end namespace ug
1086 :
1087 : #endif /* __H__UG__LIB_DISC__COMMON__GEOMETRY_UTIL__ */
|
