1 /*--------------------------------------------------------------------------*/
2 /* ALBERTA: an Adaptive multi Level finite element toolbox using */
3 /* Bisectioning refinement and Error control by Residual */
4 /* Techniques for scientific Applications */
5 /* */
6 /* file: element_2d.c */
7 /* */
8 /* */
9 /* description: routines on elements that depend on the dimension in 2d */
10 /* */
11 /*--------------------------------------------------------------------------*/
12 /* */
13 /* authors: Alfred Schmidt */
14 /* Zentrum fuer Technomathematik */
15 /* Fachbereich 3 Mathematik/Informatik */
16 /* Universitaet Bremen */
17 /* Bibliothekstr. 2 */
18 /* D-28359 Bremen, Germany */
19 /* */
20 /* Kunibert G. Siebert */
21 /* Institut fuer Mathematik */
22 /* Universitaet Augsburg */
23 /* Universitaetsstr. 14 */
24 /* D-86159 Augsburg, Germany */
25 /* */
26 /* Claus-Justus Heine */
27 /* Abteilung f�r Angewandte Mathematik */
28 /* Universitaet Freiburg */
29 /* Hermann-Herder-Str. 10 */
30 /* 79104 Freiburg, Germany */
31 /* */
32 /* http://www.alberta-fem.de */
33 /* */
34 /* (c) by A. Schmidt and K.G. Siebert (1996-2003), */
35 /* C.-J. Heine (2006-2007). */
36 /* */
37 /*--------------------------------------------------------------------------*/
38
39 #ifdef HAVE_CONFIG_H
40 # include "config.h" /* probably a good idea to pull this one in here ... */
41 #endif
42
43 #include "alberta.h"
44
45 /*--------------------------------------------------------------------------*/
46 /* world_to_coord_2d(): return -1 if inside, otherwise index of lambda<0 */
47 /*--------------------------------------------------------------------------*/
48
world_to_coord_2d(const EL_INFO * el_info,const REAL * xy,REAL_B lambda)49 int world_to_coord_2d(const EL_INFO *el_info, const REAL *xy,
50 REAL_B lambda)
51 {
52 FUNCNAME("world_to_coord_2d");
53 REAL edge[2][DIM_OF_WORLD], x[DIM_OF_WORLD];
54 REAL x0, det, det0, det1, adet, lmin;
55 int i, j, k;
56
57 #if DIM_OF_WORLD != 2
58 ERROR_EXIT("Not yet implemented for DIM_OF_WORLD != 2\n");
59 #endif
60
61 DEBUG_TEST_EXIT((el_info->fill_flag & FILL_COORDS) ||
62 !el_info->mesh->parametric ||
63 el_info->mesh->parametric->use_reference_mesh,
64 "You must enable the use_reference_mesh entry in the "
65 "PARAMETRIC structure to use this function on the "
66 "reference mesh. Use parametric->coord_to_world() "
67 "to access the parametric mesh\n");
68
69 /* wir haben das gleichungssystem zu loesen: */
70 /* ( q1x q2x ) (lambda1) (qx) */
71 /* ( q1y q2y ) (lambda2) = (qy) */
72 /* mit qi=pi-p3, q=xy-p3 */
73
74 for (j=0; j<DIM_OF_WORLD; j++) {
75 x0 = el_info->coord[2][j];
76 x[j] = xy[j] - x0;
77 for (i=0; i < 2; i++)
78 edge[i][j] = el_info->coord[i][j] - x0;
79 }
80
81 det = edge[0][0] * edge[1][1] - edge[0][1] * edge[1][0];
82 det0 = x[0] * edge[1][1] - x[1] * edge[1][0];
83 det1 = edge[0][0] * x[1] - edge[0][1] * x[0];
84
85 adet = ABS(det);
86
87 if (adet < 1.E-20) {
88 ERROR_EXIT("det = %le; abort\n", det);
89 return(-2);
90 }
91
92 lambda[0] = det0 / det;
93 lambda[1] = det1 / det;
94 lambda[2] = 1.0 - lambda[0] - lambda[1];
95
96 k = -1;
97 lmin = 0.0;
98 j = 0;
99 for (i = 0; i <= 2; i++) {
100 if ((lambda[i]*adet) < -10*REAL_EPSILON) {
101 if (lambda[i] < lmin) {
102 k = i;
103 lmin = lambda[i];
104 }
105 j++;
106 }
107 }
108
109 return(k);
110 }
111
112 /*--------------------------------------------------------------------------*/
113 /* transform local coordinates l to world coordinates; if w is non NULL */
114 /* store them at w otherwise return a pointer to some local static */
115 /* area containing world coordintes */
116 /*--------------------------------------------------------------------------*/
117
coord_to_world_2d(const EL_INFO * el_info,const REAL * l,REAL_D w)118 const REAL *coord_to_world_2d(const EL_INFO *el_info, const REAL *l, REAL_D w)
119 {
120 FUNCNAME("coord_to_world_2d");
121 static REAL_D world;
122 REAL *ret;
123 int i;
124
125 #if DIM_OF_WORLD < 2
126 # error Does not make sense for DIM_OF_WORLD < 2!
127 #endif
128
129 DEBUG_TEST_EXIT((el_info->fill_flag & FILL_COORDS) ||
130 !el_info->mesh->parametric ||
131 el_info->mesh->parametric->use_reference_mesh,
132 "You must enable the use_reference_mesh entry in the "
133 "PARAMETRIC structure to use this function on the "
134 "reference mesh. Use parametric->coord_to_world() "
135 "to access the parametric mesh\n");
136
137 ret = w ? w : world;
138
139 for(i = 0; i < DIM_OF_WORLD; i++)
140 ret[i] = el_info->coord[0][i] * l[0]
141 + el_info->coord[1][i] * l[1]
142 + el_info->coord[2][i] * l[2];
143
144
145 return((const REAL *) ret);
146 }
147
148 /*--------------------------------------------------------------------------*/
149 /* compute volume of an element */
150 /*--------------------------------------------------------------------------*/
151
el_det_2d(const EL_INFO * el_info)152 REAL el_det_2d(const EL_INFO *el_info)
153 {
154 FUNCNAME("el_det_2d");
155 REAL_D e1, e2;
156 REAL det;
157 const REAL *v0;
158 int i;
159 #if DIM_OF_WORLD == 3
160 REAL_D n;
161 #endif
162
163 DEBUG_TEST_EXIT((el_info->fill_flag & FILL_COORDS) ||
164 !el_info->mesh->parametric ||
165 el_info->mesh->parametric->use_reference_mesh,
166 "You must enable the use_reference_mesh entry in the "
167 "PARAMETRIC structure to use this function on the "
168 "reference mesh. Use parametric->coord_to_world() "
169 "to access the parametric mesh\n");
170
171 #if DIM_OF_WORLD < 2
172 # error Does not make sense for DIM_OF_WORLD < 2!
173 #endif
174
175 v0 = el_info->coord[0];
176 for (i = 0; i < DIM_OF_WORLD; i++) {
177 e1[i] = el_info->coord[1][i] - v0[i];
178 e2[i] = el_info->coord[2][i] - v0[i];
179 }
180
181 #if DIM_OF_WORLD==2
182 det = WEDGE_DOW(e1, e2);
183 det = ABS(det);
184 #elif DIM_OF_WORLD == 3
185 WEDGE_DOW(e1, e2, n);
186 det = NORM_DOW(n);
187 #else
188 /* In general: just the square-root of the det of the metric tensor. */
189 det = NRM2_DOW(e1)*NRM2_DOW(e2) - SQR(SCP_DOW(e1,e2));
190 det = sqrt(det);
191 #endif
192
193 return det;
194 }
195
el_volume_2d(const EL_INFO * el_info)196 REAL el_volume_2d(const EL_INFO *el_info)
197 {
198 return 0.5 * el_det_2d(el_info);
199 }
200
201 /*--------------------------------------------------------------------------*/
202 /* compute gradients of basis functions on element; return the absolute */
203 /* value of the determinant from the transformation to the reference */
204 /* element */
205 /* */
206 /* Notice: for dim < DIM_OF_WORLD grd_lam will contain tangential */
207 /* derivatives of the barycentric coordinates! */
208 /*--------------------------------------------------------------------------*/
209
el_grd_lambda_2d(const EL_INFO * el_info,REAL_BD grd_lam)210 REAL el_grd_lambda_2d(const EL_INFO *el_info, REAL_BD grd_lam)
211 {
212 FUNCNAME("el_grd_lambda_2d");
213 int i, j;
214 REAL_D e1, e2;
215 REAL det, adet;
216 const REAL *v0;
217 #if DIM_OF_WORLD <= 3
218 REAL a11, a12, a21, a22;
219 #endif
220 #if DIM_OF_WORLD == 3
221 REAL_D n;
222 REAL a13, a23;
223 #elif DIM_OF_WORLD > 3
224 REAL nrm2e1, nrm2e2, scpe12;
225 #endif
226
227 #if DIM_OF_WORLD < 2
228 # error Does not make sense for DIM_OF_WORLD < 2!
229 #endif
230
231 DEBUG_TEST_EXIT((el_info->fill_flag & FILL_COORDS) ||
232 !el_info->mesh->parametric ||
233 el_info->mesh->parametric->use_reference_mesh,
234 "You must enable the use_reference_mesh entry in the "
235 "PARAMETRIC structure to use this function on the "
236 "reference mesh. Use parametric->coord_to_world() "
237 "to access the parametric mesh\n");
238
239 v0 = el_info->coord[0];
240 for (i = 0; i < DIM_OF_WORLD; i++) {
241 e1[i] = el_info->coord[1][i] - v0[i];
242 e2[i] = el_info->coord[2][i] - v0[i];
243 }
244
245 #if DIM_OF_WORLD == 2
246 det = WEDGE_DOW(e1, e2);
247 adet = ABS(det);
248 #elif DIM_OF_WORLD == 3
249 WEDGE_DOW(e1, e2, n);
250 det = NRM2_DOW(n);
251 adet = sqrt(det);
252 #else
253 nrm2e1 = NRM2_DOW(e1);
254 nrm2e2 = NRM2_DOW(e2);
255 scpe12 = SCP_DOW(e1,e2);
256 det = nrm2e1*nrm2e2 - SQR(scpe12);
257 adet = sqrt(det);
258 #endif
259
260 if (adet < 1.0E-25) {
261 /* cH: is this _really_ a good test? Shouldn't the test relate to
262 * the reference triangle (i.e. shouldn't we just catch degenerate
263 * triangles)
264 */
265 MSG("abs(det) = %lf\n", adet);
266
267 for (i = 0; i < N_LAMBDA_MAX; i++)
268 for (j = 0; j < DIM_OF_WORLD; j++)
269 grd_lam[i][j] = 0.0;
270 } else {
271 det = 1.0 / det;
272 #if DIM_OF_WORLD == 2
273 a11 = e2[1] * det;
274 a21 = -e2[0] * det;
275 a12 = -e1[1] * det;
276 a22 = e1[0] * det;
277
278 grd_lam[1][0] = a11;
279 grd_lam[1][1] = a21;
280 grd_lam[2][0] = a12;
281 grd_lam[2][1] = a22;
282 grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
283 grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
284 #elif DIM_OF_WORLD == 3
285 a11 = (e2[1]* n[2] - e2[2]* n[1]) * det;
286 a12 = (e2[2]* n[0] - e2[0]* n[2]) * det;
287 a13 = (e2[0]* n[1] - e2[1]* n[0]) * det;
288 a21 = (e1[2]* n[1] - e1[1]* n[2]) * det;
289 a22 = (e1[0]* n[2] - e1[2]* n[0]) * det;
290 a23 = (e1[1]* n[0] - e1[0]* n[1]) * det;
291
292 grd_lam[1][0] = a11;
293 grd_lam[1][1] = a12;
294 grd_lam[1][2] = a13;
295 grd_lam[2][0] = a21;
296 grd_lam[2][1] = a22;
297 grd_lam[2][2] = a23;
298 grd_lam[0][0] = -grd_lam[1][0] - grd_lam[2][0];
299 grd_lam[0][1] = -grd_lam[1][1] - grd_lam[2][1];
300 grd_lam[0][2] = -grd_lam[1][2] - grd_lam[2][2];
301 #else
302 /* So what is \nabla_x \lambda_i? It is the vector orthogonal to
303 * edge number i (pointing to vertex i) with length 1/|h_i|, where
304 * h_i is the height vector between egde i and vertex i. This is
305 * most easily computed as:
306 */
307 AXPBY_DOW(nrm2e2, e1, -scpe12, e2, grd_lam[1]);
308 SCAL_DOW(det, grd_lam[1]);
309 AXPBY_DOW(nrm2e1, e2, -scpe12, e1, grd_lam[2]);
310 SCAL_DOW(det, grd_lam[2]);
311 AXPBY_DOW(-1.0, grd_lam[1], -1.0, grd_lam[2], grd_lam[0]);
312 #endif
313 }
314
315 /* Note: we have to clear any excess entries to 0, otherwise
316 * eval_grd_uh_fast() and friends will not work.
317 */
318 for (i = N_LAMBDA_2D; i < N_LAMBDA_MAX; i++) {
319 SET_DOW(0.0, grd_lam[i]);
320 }
321
322 return adet;
323 }
324
325
326 /*--------------------------------------------------------------------------*/
327 /* calculate a normal of an edge of a triangle with coordinates */
328 /* coord; return the absolute value of the determinant from the */
329 /* transformation to the reference element */
330 /*--------------------------------------------------------------------------*/
331
get_wall_normal_2d(const EL_INFO * el_info,int i0,REAL * normal)332 REAL get_wall_normal_2d(const EL_INFO *el_info, int i0, REAL *normal)
333 {
334 #if ALBERTA_DEBUG || DIM_OF_WORLD == 2
335 FUNCNAME("get_face_normal_2d");
336 #endif
337 static const int ind[5] = {0, 1, 2, 0, 1};
338 REAL det;
339 int i1 = ind[i0+1], i2 = ind[i0+2];
340 const REAL_D *coord = el_info->coord;
341 REAL_D tmp_normal, e2;
342 #if DIM_OF_WORLD > 2
343 REAL_D w;
344 #endif
345
346 DEBUG_TEST_EXIT((el_info->fill_flag & FILL_COORDS) ||
347 !el_info->mesh->parametric ||
348 el_info->mesh->parametric->use_reference_mesh,
349 "You must enable the use_reference_mesh entry in the "
350 "PARAMETRIC structure to use this function on the "
351 "reference mesh. Use parametric->coord_to_world() "
352 "to access the parametric mesh\n");
353
354 if (normal == NULL) {
355 normal = tmp_normal;
356 }
357
358 #if DIM_OF_WORLD == 2
359 normal[0] = coord[i2][1] - coord[i1][1];
360 normal[1] = coord[i1][0] - coord[i2][0];
361
362 det = NORM_DOW(normal);
363 TEST_EXIT(det > 1.e-30, "det = 0 on face %d\n", i0);
364
365 SCAL_DOW(1.0/det, normal);
366
367 /* The mesh is not necessarily oriented, so: */
368 AXPBY_DOW(1.0, coord[i2], -1.0, coord[i0], e2);
369 if (SCP_DOW(e2, normal) < 0.0) {
370 SCAL_DOW(-1.0, normal);
371 }
372
373 #else
374 AXPBY_DOW(1.0, coord[i2], -1.0, coord[i0], e2);
375 AXPBY_DOW(1.0, coord[i2], -1.0, coord[i1], w);
376
377 det = NRM2_DOW(w);
378 AXPBY_DOW(det, e2, -SCP_DOW(w, e2), w, normal);
379 SCAL_DOW(1.0/NORM_DOW(normal), normal);
380 det = sqrt(det);
381 #endif
382
383 return det;
384 }
385
386 /*--------------------------------------------------------------------------*/
387 /* orient the vertices of walls */
388 /* used by estimator for the jumps => same quadrature nodes from both sides!*/
389 /*--------------------------------------------------------------------------*/
390
391 const int sorted_wall_vertices_2d[N_WALLS_2D][DIM_FAC_2D][2*N_VERTICES_1D-1] = {
392 {{1, 2, 1}, {2, 1, 2}},
393 {{2, 0, 2}, {0, 2, 0}},
394 {{0, 1, 0}, {1, 0, 1}}
395 };
396
397 /* The return value "perm" is the index into sorted_wall_vertices_3d
398 * such that
399 *
400 * nv = sorted_wall_vertices_2d[oppv][perm][i]
401 *
402 * matches
403 *
404 * v = vertex_of_wall_2d[wall][i]
405 *
406 * I.e.: el->dof[v][0] == neigh->dof[nv][0] is true
407 *
408 * The benefit over wall_orientation_2d() is that we do not need to
409 * keep two permutations into account, but just one. This simplifies
410 * the definition of wall-quadratures.
411 */
wall_rel_orientation_2d(const EL * el,const EL * neigh,int wall,int ov)412 int wall_rel_orientation_2d(const EL *el, const EL *neigh, int wall, int ov)
413 {
414 const int *vow = vertex_of_wall_2d[wall];
415 const int *vow_n = vertex_of_wall_2d[ov];
416 DOF **dof = el->dof;
417 DOF **dof_n = neigh->dof;
418
419 return dof[vow[0]][0] != dof_n[vow_n[0]][0];
420 }
421
wall_orientation_2d(const EL * el,int wall)422 int wall_orientation_2d(const EL *el, int wall)
423 {
424 int no;
425 const int *vow = vertex_of_wall_2d[wall];
426 DOF **dof = el->dof;
427
428 no = dof[vow[0]][0] > dof[vow[1]][0];
429
430 return no;
431 }
432
433