1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5
6 #include <math.h>
7 #include "GetDPConfig.h"
8 #include "ProData.h"
9 #include "ProDefine.h"
10 #include "F.h"
11 #include "NumericUtils.h"
12 #include "MallocUtils.h"
13 #include "Message.h"
14
15 #if defined(HAVE_KERNEL)
16 #include "GeoData.h"
17 #include "DofData.h"
18 #include "Get_Geometry.h"
19 #endif
20
21 #define SQU(a) ((a)*(a))
22
23 extern struct CurrentData Current ;
24
F_Normal(F_ARG)25 void F_Normal(F_ARG)
26 {
27 #if !defined(HAVE_KERNEL)
28 Message::Error("F_Normal requires Kernel");
29 #else
30 int k ;
31
32 if(!Current.Element || Current.Element->Num == NO_ELEMENT)
33 Message::Error("No element on which to compute 'F_Normal'");
34
35 Geo_CreateNormal(Current.Element->Type,
36 Current.Element->x,
37 Current.Element->y,
38 Current.Element->z,
39 V->Val);
40
41 if (Current.NbrHar != 1) {
42 V->Val[MAX_DIM] = 0. ;
43 V->Val[MAX_DIM+1] = 0. ;
44 V->Val[MAX_DIM+2] = 0. ;
45 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
46 V->Val[MAX_DIM* k ] = V->Val[0] ;
47 V->Val[MAX_DIM* k +1] = V->Val[1] ;
48 V->Val[MAX_DIM* k +2] = V->Val[2] ;
49 V->Val[MAX_DIM*(k+1) ] = 0. ;
50 V->Val[MAX_DIM*(k+1)+1] = 0. ;
51 V->Val[MAX_DIM*(k+1)+2] = 0. ;
52 }
53 }
54 V->Type = VECTOR ;
55 #endif
56 }
57
F_NormalSource(F_ARG)58 void F_NormalSource(F_ARG)
59 {
60 #if !defined(HAVE_KERNEL)
61 Message::Error("F_NormalSource requires Kernel");
62 #else
63 int k ;
64
65 if(!Current.ElementSource || Current.ElementSource->Num == NO_ELEMENT)
66 Message::Error("No element on which to compute 'F_NormalSource'");
67
68 Geo_CreateNormal(Current.ElementSource->Type,
69 Current.ElementSource->x,
70 Current.ElementSource->y,
71 Current.ElementSource->z,
72 V->Val);
73
74 if (Current.NbrHar != 1) {
75 V->Val[MAX_DIM] = 0. ;
76 V->Val[MAX_DIM+1] = 0. ;
77 V->Val[MAX_DIM+2] = 0. ;
78 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
79 V->Val[MAX_DIM* k ] = V->Val[0] ;
80 V->Val[MAX_DIM* k +1] = V->Val[1] ;
81 V->Val[MAX_DIM* k +2] = V->Val[2] ;
82 V->Val[MAX_DIM*(k+1) ] = 0. ;
83 V->Val[MAX_DIM*(k+1)+1] = 0. ;
84 V->Val[MAX_DIM*(k+1)+2] = 0. ;
85 }
86 }
87 V->Type = VECTOR ;
88 #endif
89 }
90
F_Tangent(F_ARG)91 void F_Tangent(F_ARG)
92 {
93 int k ;
94 double tx, ty, tz, norm ;
95
96 if(!Current.Element || Current.Element->Num == NO_ELEMENT)
97 Message::Error("No element on which to compute 'F_Tangent'");
98
99 switch (Current.Element->Type) {
100
101 case LINE :
102 tx = Current.Element->x[1] - Current.Element->x[0] ;
103 ty = Current.Element->y[1] - Current.Element->y[0] ;
104 tz = Current.Element->z[1] - Current.Element->z[0] ;
105 norm = sqrt(SQU(tx)+SQU(ty)+SQU(tz)) ;
106 V->Val[0] = tx/norm ;
107 V->Val[1] = ty/norm ;
108 V->Val[2] = tz/norm ;
109 break ;
110
111 default :
112 Message::Error("Function 'Tangent' only valid for Line Elements");
113 }
114
115 if (Current.NbrHar != 1) {
116 V->Val[MAX_DIM] = 0. ;
117 V->Val[MAX_DIM+1] = 0. ;
118 V->Val[MAX_DIM+2] = 0. ;
119 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
120 V->Val[MAX_DIM* k ] = V->Val[0] ;
121 V->Val[MAX_DIM* k +1] = V->Val[1] ;
122 V->Val[MAX_DIM* k +2] = V->Val[2] ;
123 V->Val[MAX_DIM*(k+1) ] = 0. ;
124 V->Val[MAX_DIM*(k+1)+1] = 0. ;
125 V->Val[MAX_DIM*(k+1)+2] = 0. ;
126 }
127 }
128 V->Type = VECTOR ;
129 }
130
F_TangentSource(F_ARG)131 void F_TangentSource(F_ARG)
132 {
133 int k ;
134 double tx, ty, tz, norm ;
135
136 if(!Current.ElementSource || Current.ElementSource->Num == NO_ELEMENT)
137 Message::Error("No element on which to compute 'F_TangentSource'");
138
139 switch (Current.ElementSource->Type) {
140
141 case LINE :
142 tx = Current.ElementSource->x[1] - Current.ElementSource->x[0] ;
143 ty = Current.ElementSource->y[1] - Current.ElementSource->y[0] ;
144 tz = Current.ElementSource->z[1] - Current.ElementSource->z[0] ;
145 norm = sqrt(SQU(tx)+SQU(ty)+SQU(tz)) ;
146 V->Val[0] = tx/norm ;
147 V->Val[1] = ty/norm ;
148 V->Val[2] = tz/norm ;
149 break ;
150
151 default :
152 Message::Error("Function 'TangentSource' only valid for Line Elements");
153 }
154
155 if (Current.NbrHar != 1) {
156 V->Val[MAX_DIM] = 0. ;
157 V->Val[MAX_DIM+1] = 0. ;
158 V->Val[MAX_DIM+2] = 0. ;
159 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
160 V->Val[MAX_DIM* k ] = V->Val[0] ;
161 V->Val[MAX_DIM* k +1] = V->Val[1] ;
162 V->Val[MAX_DIM* k +2] = V->Val[2] ;
163 V->Val[MAX_DIM*(k+1) ] = 0. ;
164 V->Val[MAX_DIM*(k+1)+1] = 0. ;
165 V->Val[MAX_DIM*(k+1)+2] = 0. ;
166 }
167 }
168 V->Type = VECTOR ;
169 }
170
F_ElementVol(F_ARG)171 void F_ElementVol(F_ARG)
172 {
173 #if !defined(HAVE_KERNEL)
174 Message::Error("F_ElementVol requires Kernel");
175 #else
176 int k;
177 double Vol = 0.;
178 MATRIX3x3 Jac;
179
180 if(!Current.Element || Current.Element->Num == NO_ELEMENT)
181 Message::Error("No element on which to compute 'F_ElementVol'");
182
183 /* It would be more efficient to compute the volumes directly from
184 the node coordinates, but I'm lazy... */
185
186 Get_NodesCoordinatesOfElement(Current.Element) ;
187 Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w) ;
188
189 /* we use the most general case (3D embedding) */
190
191 switch(Current.Element->Type){
192 case LINE:
193 Vol = 2. * JacobianLin3D(Current.Element,&Jac);
194 break;
195 case TRIANGLE:
196 Vol = 0.5 * JacobianSur3D(Current.Element,&Jac) ;
197 break;
198 case QUADRANGLE:
199 Vol = 4. * JacobianSur3D(Current.Element,&Jac) ;
200 break;
201 case TETRAHEDRON:
202 Vol = 1./6. * JacobianVol3D(Current.Element,&Jac) ;
203 break;
204 case HEXAHEDRON:
205 Vol = 8. * JacobianVol3D(Current.Element,&Jac) ;
206 break;
207 case PRISM:
208 Vol = JacobianVol3D(Current.Element,&Jac) ;
209 break;
210 case PYRAMID:
211 Vol = 4./3. * JacobianVol3D(Current.Element,&Jac) ;
212 break;
213 default :
214 Message::Error("F_ElementVol not implemented for %s",
215 Get_StringForDefine(Element_Type, Current.Element->Type));
216 }
217
218 V->Type = SCALAR ;
219 V->Val[0] = fabs(Vol);
220 V->Val[MAX_DIM] = 0.;
221
222 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
223 V->Val[MAX_DIM* k] = V->Val[0] ;
224 V->Val[MAX_DIM* (k+1)] = 0. ;
225 }
226 #endif
227 }
228
F_SurfaceArea(F_ARG)229 void F_SurfaceArea(F_ARG)
230 {
231 #if !defined(HAVE_KERNEL)
232 Message::Error("F_SurfaceArea requires Kernel");
233 #else
234 struct Element Element ;
235
236 // check if the cache can be reused
237 if (Fct->Active) {
238 bool recompute = false;
239 if(Fct->NbrParameters != List_Nbr(Fct->Active->Case.SurfaceArea.RegionList)) {
240 recompute = true;
241 }
242 else if(Fct->NbrParameters >= 1) {
243 for(int i = 0; i < Fct->NbrParameters; i++) {
244 int num ; List_Read(Fct->Active->Case.SurfaceArea.RegionList, i, &num);
245 if(num != (int)(Fct->Para[i])) {
246 recompute = true;
247 break;
248 }
249 }
250 }
251 else if(Current.Region != Fct->Active->Case.SurfaceArea.RegionCurrent) {
252 recompute = true;
253 }
254 if(recompute) {
255 List_Delete(Fct->Active->Case.SurfaceArea.RegionList);
256 Free(Fct->Active);
257 Fct->Active = NULL;
258 }
259 }
260
261 if (!Fct->Active) {
262 Fct->Active = (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
263
264 List_T *InitialList_L = NULL;
265 if (Fct->NbrParameters >= 1) {
266 InitialList_L = List_Create(Fct->NbrParameters,1,sizeof(int));
267 for (int i = 0; i < Fct->NbrParameters; i++) {
268 int Index_Region = (int)(Fct->Para[i]) ;
269 List_Add(InitialList_L, &Index_Region);
270 }
271 }
272
273 double Val_Surface = 0. ;
274 int Nbr_Element = Geo_GetNbrGeoElements() ;
275 int Nbr_Found_Element = 0;
276 for (int i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {
277 Element.GeoElement = Geo_GetGeoElement(i_Element) ;
278 if ((InitialList_L &&
279 List_Search(InitialList_L, &(Element.GeoElement->Region), fcmp_int)) ||
280 (!InitialList_L && Element.GeoElement->Region == Current.Region)) {
281
282 Nbr_Found_Element++;
283
284 Element.Num = Element.GeoElement->Num ;
285 Element.Type = Element.GeoElement->Type ;
286
287 if (Element.Type == TRIANGLE || Element.Type == QUADRANGLE ||
288 Element.Type == TRIANGLE_2) {
289
290 Get_NodesCoordinatesOfElement(&Element) ;
291 Get_BFGeoElement(&Element, 0., 0., 0.) ;
292 double c11 = 0., c21 = 0., c12 = 0., c22 = 0., c13 = 0., c23 = 0.;
293 for (int i = 0 ; i < Element.GeoElement->NbrNodes ; i++ ) {
294 c11 += Element.x[i] * Element.dndu[i][0] ;
295 c21 += Element.x[i] * Element.dndu[i][1] ;
296 c12 += Element.y[i] * Element.dndu[i][0] ;
297 c22 += Element.y[i] * Element.dndu[i][1] ;
298 c13 += Element.z[i] * Element.dndu[i][0] ;
299 c23 += Element.z[i] * Element.dndu[i][1] ;
300 }
301 double DetJac = sqrt(SQU(c11 * c22 - c12 * c21) +
302 SQU(c13 * c21 - c11 * c23) +
303 SQU(c12 * c23 - c13 * c22));
304
305 if (Element.Type == TRIANGLE || Element.Type == TRIANGLE_2)
306 Val_Surface += fabs(DetJac) * 0.5 ;
307 else if (Element.Type == QUADRANGLE)
308 Val_Surface += fabs(DetJac) * 4. ;
309
310 }
311 else if (Element.Type == LINE || Element.Type == LINE_2) {
312 Get_NodesCoordinatesOfElement(&Element) ;
313 Get_BFGeoElement(&Element, 0., 0., 0.) ;
314
315 double c11 = 0., c12 = 0., c13 = 0.;
316 for (int i = 0; i < Element.GeoElement->NbrNodes; i++) {
317 c11 += Element.x[i] * Element.dndu[i][0] ;
318 c12 += Element.y[i] * Element.dndu[i][0] ;
319 c13 += Element.z[i] * Element.dndu[i][0] ;
320 }
321 double DetJac = sqrt(SQU(c11)+SQU(c12)+SQU(c13));
322
323 Val_Surface += fabs(DetJac) * 2 ; // SurfaceArea of LINE x 1m
324 }
325 else {
326 Message::Error("Function 'SurfaceArea' only valid for line, triangle or "
327 "quandrangle elements");
328 }
329 }
330 }
331 Fct->Active->Case.SurfaceArea.RegionList = InitialList_L ;
332 Fct->Active->Case.SurfaceArea.RegionCurrent = Current.Region ;
333 Fct->Active->Case.SurfaceArea.Value = Val_Surface ;
334
335 if(!Nbr_Found_Element)
336 Message::Info("No element found in SurfaceArea[]");
337 }
338
339 V->Type = SCALAR ;
340 V->Val[0] = Fct->Active->Case.SurfaceArea.Value ;
341 V->Val[MAX_DIM] = 0.;
342
343 for (int k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2) {
344 V->Val[MAX_DIM* k] = V->Val[0] ;
345 V->Val[MAX_DIM* (k+1)] = 0. ;
346 }
347 #endif
348 }
349
F_GetVolume(F_ARG)350 void F_GetVolume(F_ARG)
351 {
352 #if !defined(HAVE_KERNEL)
353 Message::Error("F_GetVolume requires Kernel");
354 #else
355 struct Element Element ;
356 List_T * InitialList_L;
357
358 int Nbr_Element, i_Element ;
359 double Val_Volume ;
360 double c11, c21, c31, c12, c22, c32, c13, c23, c33 ;
361 double DetJac ;
362 int i, k ;
363
364 if (!Fct->Active) {
365 Fct->Active = (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
366
367 if (Fct->NbrParameters == 1) {
368 int Index_Region = (int)(Fct->Para[0]) ;
369 InitialList_L = List_Create(1, 1, sizeof(int));
370 List_Add(InitialList_L,&Index_Region);
371 }
372 else {
373 InitialList_L = NULL ;
374 }
375
376 Val_Volume = 0. ;
377 Nbr_Element = Geo_GetNbrGeoElements() ;
378 for (i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {
379 Element.GeoElement = Geo_GetGeoElement(i_Element) ;
380 if ((InitialList_L &&
381 List_Search(InitialList_L, &(Element.GeoElement->Region), fcmp_int)) ||
382 (!InitialList_L && Element.GeoElement->Region == Current.Region)) {
383 Element.Num = Element.GeoElement->Num ;
384 Element.Type = Element.GeoElement->Type ;
385
386 if (Element.Type == TETRAHEDRON || Element.Type == TETRAHEDRON_2 ||
387 Element.Type == HEXAHEDRON ||
388 Element.Type == PRISM) {
389
390 Get_NodesCoordinatesOfElement(&Element) ;
391 Get_BFGeoElement(&Element, 0., 0., 0.) ;
392
393 c11 = c21 = c31 = c12 = c22 = c32 = c13 = c23 = c33 = 0;
394 for ( i = 0 ; i < Element.GeoElement->NbrNodes ; i++ ) {
395 c11 += Element.x[i] * Element.dndu[i][0] ;
396 c21 += Element.x[i] * Element.dndu[i][1] ;
397 c31 += Element.x[i] * Element.dndu[i][2] ;
398
399 c12 += Element.y[i] * Element.dndu[i][0] ;
400 c22 += Element.y[i] * Element.dndu[i][1] ;
401 c32 += Element.y[i] * Element.dndu[i][2] ;
402
403 c13 += Element.z[i] * Element.dndu[i][0] ;
404 c23 += Element.z[i] * Element.dndu[i][1] ;
405 c33 += Element.z[i] * Element.dndu[i][2] ;
406 }
407
408 DetJac = c11 * c22 * c33 + c13 * c21 * c32
409 + c12 * c23 * c31 - c13 * c22 * c31
410 - c11 * c23 * c32 - c12 * c21 * c33 ;
411
412 switch(Element.Type){
413 case TETRAHEDRON:
414 case TETRAHEDRON_2:
415 Val_Volume += 1./6. * fabs(DetJac);
416 break;
417 case HEXAHEDRON:
418 Val_Volume += 8. * fabs(DetJac);
419 break;
420 case PRISM:
421 Val_Volume += fabs(DetJac);
422 break;
423 }
424 }
425 else {
426 Message::Error("Function 'GetVolume' not valid for %s",
427 Get_StringForDefine(Element_Type, Element.Type));
428 }
429 }
430 }
431 Fct->Active->Case.GetVolume.Value = Val_Volume ;
432 }
433
434 V->Type = SCALAR ;
435 V->Val[0] = Fct->Active->Case.GetVolume.Value ;
436 V->Val[MAX_DIM] = 0.;
437
438 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
439 V->Val[MAX_DIM* k] = V->Val[0] ;
440 V->Val[MAX_DIM* (k+1)] = 0. ;
441 }
442 #endif
443 }
444
F_GetNumElement(F_ARG)445 void F_GetNumElement(F_ARG)
446 {
447 int k;
448
449 V->Type = SCALAR ;
450 V->Val[0] = (double)Current.Element->Num;
451 V->Val[MAX_DIM] = 0.;
452
453 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
454 V->Val[MAX_DIM* k] = V->Val[0] ;
455 V->Val[MAX_DIM* (k+1)] = 0 ;
456 }
457 }
458
F_GetNumElements(F_ARG)459 void F_GetNumElements(F_ARG)
460 {
461 #if !defined(HAVE_KERNEL)
462 Message::Error("F_GetNumElements requires Kernel");
463 #else
464 struct Element Element ;
465
466 if (!Fct->Active) {
467 Fct->Active = (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
468 List_T * InitialList_L;
469 if (Fct->NbrParameters == 1) {
470 int Index_Region = (int)(Fct->Para[0]) ;
471 InitialList_L = List_Create(1, 1, sizeof(int));
472 List_Add(InitialList_L, &Index_Region);
473 }
474 else if (Fct->NbrParameters > 1) {
475 InitialList_L = List_Create(Fct->NbrParameters,1,sizeof(int));
476 List_Reset(InitialList_L);
477 for (int i=0; i<Fct->NbrParameters; i++) {
478 int Index_Region = (int)(Fct->Para[i]) ;
479 List_Add(InitialList_L, &Index_Region);
480 }
481 }
482 else {
483 InitialList_L = NULL ;
484 }
485
486 int Count = 0. ;
487 int Nbr_Element = Geo_GetNbrGeoElements() ;
488 if(!InitialList_L){
489 Count = Nbr_Element ;
490 }
491 else{
492 for (int i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {
493 Element.GeoElement = Geo_GetGeoElement(i_Element) ;
494 if ((InitialList_L &&
495 List_Search(InitialList_L, &(Element.GeoElement->Region), fcmp_int)) ||
496 (!InitialList_L && Element.GeoElement->Region == Current.Region)) {
497 Count++;
498 }
499 }
500 }
501 Fct->Active->Case.GetNumElements.Value = Count ;
502 }
503
504 V->Type = SCALAR ;
505 V->Val[0] = Fct->Active->Case.GetNumElements.Value ;
506 V->Val[MAX_DIM] = 0.;
507
508 for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
509 V->Val[MAX_DIM* k] = V->Val[0] ;
510 V->Val[MAX_DIM* (k+1)] = 0. ;
511 }
512 #endif
513 }
514
F_GetNumNodes(F_ARG)515 void F_GetNumNodes(F_ARG)
516 {
517 #if !defined(HAVE_KERNEL)
518 Message::Error("F_GetNumNodes requires Kernel");
519 #else
520 // TODO: accept arguments to limit to some regions
521 V->Type = SCALAR ;
522 V->Val[0] = Geo_GetNbrGeoNodes() ;
523 V->Val[MAX_DIM] = 0.;
524 for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
525 V->Val[MAX_DIM* k] = V->Val[0] ;
526 V->Val[MAX_DIM* (k+1)] = 0. ;
527 }
528 #endif
529 }
530
F_CellSize(F_ARG)531 void F_CellSize(F_ARG)
532 {
533 #if !defined(HAVE_KERNEL)
534 Message::Error("F_CellSize requires Kernel");
535 #else
536 double cellSize, Vol;
537 MATRIX3x3 Jac;
538 double c11, c21, c12, c22, DetJac;
539 int i, k ;
540
541 if(!Current.Element || Current.Element->Num == NO_ELEMENT)
542 Message::Error("No element on which to compute 'CellSize'");
543
544 Get_NodesCoordinatesOfElement(Current.Element) ;
545 Get_BFGeoElement(Current.Element, 0., 0., 0.) ;
546
547 switch(Current.Element->Type){
548 case LINE:
549 cellSize = 2. * JacobianLin3D(Current.Element,&Jac);
550 break;
551 case TRIANGLE:
552 c11 = c21 = c12 = c22 = 0. ;
553
554 for ( i = 0 ; i < Current.Element->GeoElement->NbrNodes ; i++ ) {
555 c11 += Current.Element->x[i] * Current.Element->dndu[i][0] ;
556 c21 += Current.Element->x[i] * Current.Element->dndu[i][1] ;
557 c12 += Current.Element->y[i] * Current.Element->dndu[i][0] ;
558 c22 += Current.Element->y[i] * Current.Element->dndu[i][1] ;
559 }
560 DetJac = c11 * c22 - c12 * c21 ;
561 cellSize =
562 sqrt(SQU(Current.Element->x[1]-Current.Element->x[0])
563 +SQU(Current.Element->y[1]-Current.Element->y[0])
564 +SQU(Current.Element->z[1]-Current.Element->z[0]))
565 *
566 sqrt(SQU(Current.Element->x[2]-Current.Element->x[1])
567 +SQU(Current.Element->y[2]-Current.Element->y[1])
568 +SQU(Current.Element->z[2]-Current.Element->z[1]))
569 *
570 sqrt(SQU(Current.Element->x[0]-Current.Element->x[2])
571 +SQU(Current.Element->y[0]-Current.Element->y[2])
572 +SQU(Current.Element->z[0]-Current.Element->z[2]))
573 / fabs(DetJac) ;
574 break;
575 case QUADRANGLE:
576 // Message::Warning("Function CellSize not ready for QUADRANGLE") ;
577 Vol = 4. * JacobianSur3D(Current.Element,&Jac) ;
578 cellSize = sqrt(Vol);
579 break;
580 case TETRAHEDRON:
581 cellSize = 0.;
582 Message::Warning("Function CellSize not ready for TETRAHEDRON") ;
583 break;
584 case HEXAHEDRON:
585 cellSize = 0.;
586 Message::Warning("Function CellSize not ready for HEXAHEDRON") ;
587 break;
588 case PRISM:
589 cellSize = 0.;
590 Message::Warning("Function CellSize not ready for PRISM") ;
591 break;
592 default :
593 cellSize = 0.;
594 Message::Error("Function 'CellSize' not valid for %s",
595 Get_StringForDefine(Element_Type, Current.Element->Type));
596 }
597
598 V->Type = SCALAR ;
599 V->Val[0] = cellSize ;
600 V->Val[MAX_DIM] = 0. ;
601 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
602 V->Val[MAX_DIM* k] = V->Val[0] ;
603 V->Val[MAX_DIM* (k+1)] = 0. ;
604 }
605 #endif
606 }
607
F_SquNormEdgeValues(F_ARG)608 void F_SquNormEdgeValues(F_ARG)
609 {
610 #if !defined(HAVE_KERNEL)
611 Message::Error("F_SquNormEdgeValues requires Kernel");
612 #else
613 struct Geo_Element *GeoElement;
614 int i, *NumNodes;
615 double x [NBR_MAX_NODES_IN_ELEMENT] ;
616 double y [NBR_MAX_NODES_IN_ELEMENT] ;
617 double z [NBR_MAX_NODES_IN_ELEMENT] ;
618 double xe[2], ye[2], ze[2];
619
620 int numDofData, Code_BasisFunction, CodeExist = 0;
621 struct Dof * Dof_P = NULL;
622 double Val_Dof, Val_Dof_i, valSum, sizeEdge ;
623 int k;
624
625 if(!Current.Element || Current.Element->Num == NO_ELEMENT)
626 Message::Error("No element on which to compute 'SquNormEdgeValues'");
627
628 numDofData = (int)Fct->Para[0];
629 Code_BasisFunction = (int)Fct->Para[1];
630
631 GeoElement = Current.Element->GeoElement;
632 Geo_GetNodesCoordinates
633 (GeoElement->NbrNodes, GeoElement->NumNodes, x, y, z) ;
634
635 valSum = 0.;
636
637 if(!GeoElement->NbrEdges) Geo_CreateEdgesOfElement(GeoElement) ;
638 for(i=0 ; i<GeoElement->NbrEdges ; i++){
639 NumNodes = Geo_GetNodesOfEdgeInElement(GeoElement, i) ;
640 xe[0] = x[abs(NumNodes[0])-1]; xe[1] = x[abs(NumNodes[1])-1];
641 ye[0] = y[abs(NumNodes[0])-1]; ye[1] = y[abs(NumNodes[1])-1];
642 ze[0] = z[abs(NumNodes[0])-1]; ze[1] = z[abs(NumNodes[1])-1];
643
644 CodeExist =
645 ((Dof_P =
646 Dof_GetDofStruct(Current.DofData_P0+ numDofData,
647 Code_BasisFunction, abs(GeoElement->NumEdges[i]), 0))
648 != NULL) ;
649
650 if (CodeExist) {
651 sizeEdge = sqrt( SQU(xe[1]-xe[0]) + SQU(ye[1]-ye[0]) + SQU(ze[1]-ze[0]) );
652
653 if(Current.NbrHar==1){
654 Dof_GetRealDofValue(Current.DofData_P0+ numDofData, Dof_P, &Val_Dof) ;
655 Val_Dof = Val_Dof/sizeEdge;
656 valSum += SQU(Val_Dof) * sizeEdge;
657 }
658 else{
659 for (k = 0 ; k < Current.NbrHar ; k+=2) {
660 Dof_GetComplexDofValue
661 (Current.DofData_P0+ numDofData,
662 Dof_P + k/2*gCOMPLEX_INCREMENT,
663 &Val_Dof, &Val_Dof_i) ;
664 Val_Dof = Val_Dof /sizeEdge;
665 Val_Dof_i = Val_Dof_i/sizeEdge;
666 valSum += (SQU(Val_Dof)+SQU(Val_Dof_i)) * sizeEdge;
667 }
668 }
669 }
670
671 }
672
673 V->Type = SCALAR ;
674 V->Val[0] = valSum ;
675 V->Val[MAX_DIM] = 0. ;
676
677 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
678 V->Val[MAX_DIM* k] = V->Val[0] ;
679 V->Val[MAX_DIM* (k+1)] = 0. ;
680 }
681 #endif
682 }
683
684 static double POINT_TO_PROJECT[3], ELLIPSE_PARAMETERS[2];
685
dist_ellipse(double t)686 static double dist_ellipse(double t)
687 {
688 double x, y;
689 x = ELLIPSE_PARAMETERS[0] * cos(t);
690 y = ELLIPSE_PARAMETERS[1] * sin(t);
691 return sqrt(SQU(x - POINT_TO_PROJECT[0]) + SQU(y - POINT_TO_PROJECT[1]));
692 }
693
F_ProjectPointOnEllipse(F_ARG)694 void F_ProjectPointOnEllipse(F_ARG)
695 {
696 int k;
697 double t1 = 0., t2 = 1.e-6, t3, f1, f2, f3, tol = 1.e-4;
698 double t, x, y;
699
700 POINT_TO_PROJECT[0] = A->Val[0];
701 POINT_TO_PROJECT[1] = A->Val[1];
702 POINT_TO_PROJECT[2] = A->Val[2];
703
704 ELLIPSE_PARAMETERS[0] = Fct->Para[0] ;
705 ELLIPSE_PARAMETERS[1] = Fct->Para[1] ;
706
707 mnbrak(&t1, &t2, &t3, &f1, &f2, &f3, dist_ellipse);
708
709 if(t1 > t2){
710 t = t1;
711 t1 = t3;
712 t3 = t;
713 }
714
715 brent(t1, t2, t3, dist_ellipse, tol, &t);
716
717 x = ELLIPSE_PARAMETERS[0] * cos(t);
718 y = ELLIPSE_PARAMETERS[1] * sin(t);
719
720 /* printf("SL(%g,%g,0,%g,%g,0){1,1};\n", A->Val[0], A->Val[1], x, y); */
721
722 for (k = 0 ; k < Current.NbrHar ; k++) {
723 V->Val[MAX_DIM*k ] = 0. ;
724 V->Val[MAX_DIM*k+1] = 0. ;
725 V->Val[MAX_DIM*k+2] = 0. ;
726 }
727 V->Val[0] = x;
728 V->Val[1] = y;
729 V->Type = VECTOR ;
730 }
731