1 /****************************************************************************/
2 /* */
3 /* File: shapes.c */
4 /* */
5 /* Purpose: shape functions for triangles and quadrilaterals */
6 /* */
7 /* Author: Peter Bastian */
8 /* Interdisziplinaeres Zentrum fuer Wissenschaftliches Rechnen */
9 /* Universitaet Heidelberg */
10 /* Im Neuenheimer Feld 368 */
11 /* 6900 Heidelberg */
12 /* internet: ug@ica3.uni-stuttgart.de */
13 /* */
14 /* History: 08.04.92 begin, ug version 2.0 */
15 /* 20.11.94 moved shapes.c from ug to cd folder */
16 /* */
17 /* Remarks: */
18 /* */
19 /****************************************************************************/
20
21 /****************************************************************************/
22 /* */
23 /* include files */
24 /* system include files */
25 /* application include files */
26 /* */
27 /****************************************************************************/
28
29 #include <config.h>
30 #include <cmath>
31 #include <cassert>
32 #include <cstddef>
33
34 #include <dune/uggrid/low/architecture.h>
35 #include <dune/uggrid/low/misc.h>
36 #include <dune/uggrid/low/ugtypes.h>
37 #include "gm.h"
38 #include "evm.h"
39 #include "shapes.h"
40
41
42 USING_UG_NAMESPACE
43 USING_UGDIM_NAMESPACE
44
45 /****************************************************************************/
46 /* */
47 /* defines in the following order */
48 /* */
49 /* compile time constants defining static data size (i.e. arrays) */
50 /* other constants */
51 /* macros */
52 /* */
53 /****************************************************************************/
54
55 #define SMALL_DIFF 1e-20
56 #define MAX_ITER 20
57
58 /****************************************************************************/
59 /* */
60 /* data structures used in this source file (exported data structures are */
61 /* in the corresponding include file!) */
62 /* */
63 /****************************************************************************/
64
65 /****************************************************************************/
66 /* */
67 /* definition of exported global variables */
68 /* */
69 /****************************************************************************/
70
71 /****************************************************************************/
72 /* */
73 /* definition of variables global to this source file only (static!) */
74 /* */
75 /****************************************************************************/
76
77 /* local midpoints */
78 #ifdef __TWODIM__
79 static DOUBLE_VECTOR_2D LMP_Triangle = {0.333333333333333333, 0.333333333333333333};
80 static DOUBLE_VECTOR_2D LMP_Quadrilateral = {0.5, 0.5};
81 #endif
82 #ifdef __THREEDIM__
83 static DOUBLE_VECTOR_3D LMP_Tetrahedron = {0.25, 0.25, 0.25};
84 static DOUBLE_VECTOR_3D LMP_Pyramid = {0.5, 0.5, 0.33333333333333333};
85 static DOUBLE_VECTOR_3D LMP_Prism = {0.333333333333333333,
86 0.333333333333333333,0.5};
87 static DOUBLE_VECTOR_3D LMP_Hexahedron = {0.5, 0.5, 0.5};
88 #endif
89
90 /****************************************************************************/
91 /* */
92 /* forward declarations of functions used before they are defined */
93 /* */
94 /****************************************************************************/
95
96
97 /****************************************************************************/
98 /** \brief Local midpoint
99
100 \param n - number of corners of the element
101
102 This function gives the local coordinates of the midpoint of an element
103
104 \return
105 Pointer to the coordinate array
106 */
107 /****************************************************************************/
108
LMP(INT n)109 DOUBLE * NS_DIM_PREFIX LMP (INT n)
110 {
111 #ifdef __TWODIM__
112 switch (n)
113 {
114 case 3 : return (LMP_Triangle);
115 case 4 : return (LMP_Quadrilateral);
116 }
117 #endif
118
119 #ifdef __THREEDIM__
120 switch (n)
121 {
122 case 4 : return (LMP_Tetrahedron);
123 case 5 : return (LMP_Pyramid);
124 case 6 : return (LMP_Prism);
125 case 8 : return (LMP_Hexahedron);
126 }
127 #endif
128 return (NULL);
129 }
130
131 /****************************************************************************/
132 /** \brief Transform global coordinates to local
133
134 \param n - number of corners
135 \param Corners - coordinates of corners
136 \param EvalPoint - global coordinates
137 \param LocalCoord - local coordinates
138
139 This function transforms global coordinates to local in an evaluated point
140 in 3D.
141
142 \return <ul>
143 <li> 0 if ok </li>
144 <li> 1 if error occured. </li>
145 </ul>
146 */
147 /****************************************************************************/
148
UG_GlobalToLocal(INT n,const DOUBLE ** Corners,const DOUBLE * EvalPoint,DOUBLE * LocalCoord)149 INT NS_DIM_PREFIX UG_GlobalToLocal (INT n, const DOUBLE **Corners,
150 const DOUBLE *EvalPoint, DOUBLE *LocalCoord)
151 {
152 DOUBLE_VECTOR tmp,diff,M[DIM],IM[DIM];
153 DOUBLE s,IMdet;
154 INT i;
155
156 V_DIM_SUBTRACT(EvalPoint,Corners[0],diff);
157 if (n == DIM+1)
158 {
159 TRANSFORMATION(DIM+1,Corners,LocalCoord,M);
160 M_DIM_INVERT(M,IM,IMdet);
161 if (IMdet==0) return (2);
162 MT_TIMES_V_DIM(IM,diff,LocalCoord);
163 return(0);
164 }
165 V_DIM_CLEAR(LocalCoord);
166 TRANSFORMATION(n,Corners,LocalCoord,M);
167 M_DIM_INVERT(M,IM,IMdet);
168 if (IMdet==0) return (3);
169 MT_TIMES_V_DIM(IM,diff,LocalCoord);
170 for (i=0; i<MAX_ITER; i++)
171 {
172 LOCAL_TO_GLOBAL (n,Corners,LocalCoord,tmp);
173 V_DIM_SUBTRACT(tmp,EvalPoint,diff);
174 V_DIM_EUKLIDNORM(diff,s);
175 PRINTDEBUG(gm,1,("UG_GlobalToLocal %d %g\n",i,s));
176 if (s * s <= SMALL_DIFF * IMdet)
177 return (0);
178 TRANSFORMATION(n,Corners,LocalCoord,M);
179 M_DIM_INVERT(M,IM,IMdet);
180 if (IMdet==0) return (4);
181 MT_TIMES_V_DIM(IM,diff,tmp);
182 V_DIM_SUBTRACT(LocalCoord,tmp,LocalCoord);
183 }
184
185 return(1);
186 }
187
188 /****************************************************************************/
189 /** \brief Calculate inner normals of tetrahedra
190
191 \param theCorners - list of pointers to phys corner vectors
192 \param theNormals - normals of tetrahedra
193
194 This function calculates the inner normals on the sides of a tetrahedron.
195
196 \return <ul>
197 <li> 0 if ok </li>
198 <li> 1 if error occured. </li>
199 </ul>
200 */
201 /****************************************************************************/
202
203 #ifdef __THREEDIM__
TetraSideNormals(ELEMENT * theElement,DOUBLE ** theCorners,DOUBLE_VECTOR theNormals[MAX_SIDES_OF_ELEM])204 INT NS_DIM_PREFIX TetraSideNormals (ELEMENT *theElement, DOUBLE **theCorners, DOUBLE_VECTOR theNormals[MAX_SIDES_OF_ELEM])
205 {
206 ELEMENT e;
207 DOUBLE_VECTOR a, b;
208 DOUBLE h;
209 INT j,k;
210
211 /* TODO: changed MAX_CORNERS_OF_ELEM to 4 and subsequently*/
212 SETTAG(&e,4);
213 for (j=0; j<4; j++)
214 {
215 k = SIDE_OPP_TO_CORNER(&e,j);
216 V3_SUBTRACT(theCorners[(j+1)%4],theCorners[(j+2)%4],a)
217 V3_SUBTRACT(theCorners[(j+1)%4],theCorners[(j+3)%4],b)
218 V3_VECTOR_PRODUCT(a,b,theNormals[k])
219 V3_Normalize(theNormals[k]);
220 V3_SUBTRACT(theCorners[j],theCorners[(j+1)%4],a)
221 V3_SCALAR_PRODUCT(theNormals[k],a,h);
222 if (std::abs(h)<SMALL_C) return (1);
223 if (h<0.0)
224 V3_SCALE(-1.0,theNormals[k]);
225 }
226
227 return (0);
228
229 }
230 #endif
231
232 /****************************************************************************/
233 /** \brief Calculate maximal side angle of Tetrahedron
234
235 \param theCorners - list of pointers to phys corner vectors
236 \param MaxAngle - max side angle
237
238 This function calculates the maximal side angle of the tetrahedron.
239
240 \return <ul>
241 <li> 0 if ok </li>
242 <li> 1 if error occured. </li>
243 </ul>
244 */
245 /****************************************************************************/
246
247 #ifdef __THREEDIM__
TetMaxSideAngle(ELEMENT * theElement,const DOUBLE ** theCorners,DOUBLE * MaxAngle)248 INT NS_DIM_PREFIX TetMaxSideAngle (ELEMENT *theElement, const DOUBLE **theCorners, DOUBLE *MaxAngle)
249 {
250 DOUBLE_VECTOR theNormal[MAX_SIDES_OF_ELEM];
251 DOUBLE max,help;
252 INT i;
253
254 if (TetraSideNormals (theElement,(DOUBLE **)theCorners,theNormal)) return (1);
255 max = -1.0;
256 for (i=0; i<EDGES_OF_ELEM(theElement); i++)
257 {
258 V3_SCALAR_PRODUCT(theNormal[SIDE_WITH_EDGE(theElement,i,0)],theNormal[SIDE_WITH_EDGE(theElement,i,1)],help)
259 max = MAX(help,max);
260 }
261 max = MIN(max,1.0);
262 *MaxAngle = 180.0/PI*acos(-max);
263
264 return (0);
265 }
266 #endif
267
268 /****************************************************************************/
269 /** \brief Calculates side angle and length of edge of Tetrahedron
270
271 \param theCorners - list of pointers to phys corner vectors
272 \param Angle - side angle
273 \param Length - sidelength
274
275 This function calculates the side angle of a tetrahedron and the lengths of
276 the edges belonging to this side.
277
278 \return <ul>
279 <li> 0 if ok </li>
280 <li> 1 if error occured. </li>
281 </ul>
282 */
283 /****************************************************************************/
284
285 #ifdef __THREEDIM__
TetAngleAndLength(ELEMENT * theElement,const DOUBLE ** theCorners,DOUBLE * Angle,DOUBLE * Length)286 INT NS_DIM_PREFIX TetAngleAndLength (ELEMENT *theElement, const DOUBLE **theCorners, DOUBLE *Angle, DOUBLE *Length)
287 {
288 DOUBLE_VECTOR theNormals[MAX_SIDES_OF_ELEM],theEdge[MAX_EDGES_OF_ELEM];
289 DOUBLE h;
290 INT j,k;
291
292 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
293 {
294 V3_SUBTRACT(theCorners[CORNER_OF_EDGE(theElement,j,1)],theCorners[CORNER_OF_EDGE(theElement,j,0)],theEdge[j])
295 V3_EUKLIDNORM(theEdge[j],Length[j])
296 }
297 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
298 {
299 V3_VECTOR_PRODUCT(theEdge[EDGE_OF_SIDE(theElement,j,0)],theEdge[EDGE_OF_SIDE(theElement,j,1)],theNormals[j])
300 V3_Normalize(theNormals[j]);
301 k = EDGE_OF_CORNER(theElement,CORNER_OPP_TO_SIDE(theElement,j),0);
302 V3_SCALAR_PRODUCT(theNormals[j],theEdge[k],h)
303 if (std::abs(h)<SMALL_C) return (1);
304 if ( (h<0.0 && CORNER_OF_EDGE(theElement,k,1)==CORNER_OPP_TO_SIDE(theElement,j)) ||
305 (h>0.0 && CORNER_OF_EDGE(theElement,k,0)==CORNER_OPP_TO_SIDE(theElement,j)) )
306 V3_SCALE(-1.0,theNormals[j]);
307 }
308 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
309 {
310 V3_SCALAR_PRODUCT(theNormals[SIDE_WITH_EDGE(theElement,j,0)],theNormals[SIDE_WITH_EDGE(theElement,j,1)],Angle[j])
311 Angle[j] = MAX(Angle[j],-1.0);
312 Angle[j] = MIN(Angle[j], 1.0);
313 Angle[j] = (DOUBLE)acos((double)Angle[j]);
314 }
315
316 return (0);
317 }
318 #endif
319