Skip to content

Commit 355edce

Browse files
softmattertheorysoftmattertheory
authored
Merge pull request #301 from Morpho-lang/cg2-3d
CG2 elements in 3D.
2 parents 700b6c2 + e64d66b commit 355edce

6 files changed

Lines changed: 319 additions & 17 deletions

File tree

‎src/geometry/fespace.c‎

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -314,11 +314,176 @@ fespace cg2_2d = {
314314
.lower = cg2_2d_lower
315315
};
316316

317+
/* -------------------------------------------------------
318+
* CG1 element in 3D
319+
* ------------------------------------------------------- */
320+
321+
/* z=0 z=1
322+
* 2
323+
* |\
324+
* 0-1 3 // One degree of freedom per vertex
325+
*/
326+
327+
void cg1_3dinterpolate(double *lambda, double *wts) {
328+
wts[0]=lambda[0];
329+
wts[1]=lambda[1];
330+
wts[2]=lambda[2];
331+
wts[3]=lambda[3];
332+
}
333+
334+
void cg1_3dgrad(double *lambda, double *grad) {
335+
double g[] =
336+
{ 1, 0, 0, 0,
337+
0, 1, 0, 0,
338+
0, 0, 1, 0,
339+
0, 0, 0, 1 };
340+
memcpy(grad, g, sizeof(g));
341+
}
342+
343+
unsigned int cg1_3dshape[] = { 1, 0, 0, 0 };
344+
345+
double cg1_3dnodes[] = { 0.0, 0.0, 0.0,
346+
1.0, 0.0, 0.0,
347+
0.0, 1.0, 0.0,
348+
0.0, 0.0, 1.0 };
349+
350+
eldefninstruction cg1_3deldefn[] = {
351+
QUANTITY(0,0,0), // Fetch quantity on vertex 0
352+
QUANTITY(0,1,0), // Fetch quantity on vertex 1
353+
QUANTITY(0,2,0), // Fetch quantity on vertex 2
354+
QUANTITY(0,3,0), // Fetch quantity on vertex 3
355+
ENDDEFN
356+
};
357+
358+
fespace *cg1_3d_lower[] = {
359+
&cg1_2d,
360+
&cg1_1d,
361+
NULL
362+
};
363+
364+
fespace cg1_3d = {
365+
.name = "CG1",
366+
.grade = 3,
367+
.shape = cg1_3dshape,
368+
.degree = 1,
369+
.nnodes = 4,
370+
.nsubel = 0,
371+
.nodes = cg1_3dnodes,
372+
.ifn = cg1_3dinterpolate,
373+
.gfn = cg1_3dgrad,
374+
.eldefn = cg1_3deldefn,
375+
.lower = cg1_3d_lower
376+
};
377+
378+
/* -------------------------------------------------------
379+
* CG2 element in 3D
380+
* ------------------------------------------------------- */
381+
382+
/* z=0 z=0.5 z=1
383+
* 2
384+
* |\
385+
* 6 5 9
386+
* | \ | \
387+
* 0-4-1 7--8 3 - i.e. vertices
388+
*/
389+
390+
void cg2_3dinterpolate(double *lambda, double *wts) {
391+
wts[0]=lambda[0]*(2*lambda[0]-1);
392+
wts[1]=lambda[1]*(2*lambda[1]-1);
393+
wts[2]=lambda[2]*(2*lambda[2]-1);
394+
wts[3]=lambda[3]*(2*lambda[3]-1);
395+
wts[4]=4*lambda[0]*lambda[1];
396+
wts[5]=4*lambda[1]*lambda[2];
397+
wts[6]=4*lambda[2]*lambda[0];
398+
wts[7]=4*lambda[0]*lambda[3];
399+
wts[8]=4*lambda[1]*lambda[3];
400+
wts[9]=4*lambda[2]*lambda[3];
401+
}
402+
403+
void cg2_3dgrad(double *lambda, double *grad) { // TODO: FIX
404+
// Gij = d Xi[i] / d lambda[j]
405+
// Note this is in column-major order!
406+
double g[] =
407+
{ 4*lambda[0]-1, 0, 0, 0,
408+
4*lambda[1], 0, 4*lambda[2], 4*lambda[3], 0, 0,
409+
410+
0, 4*lambda[1]-1, 0, 0,
411+
4*lambda[0], 4*lambda[2], 0, 0, 4*lambda[3], 0,
412+
413+
0, 0, 4*lambda[2]-1, 0,
414+
0, 4*lambda[1], 4*lambda[0], 0, 0, 4*lambda[3],
415+
416+
0, 0, 0, 4*lambda[3]-1,
417+
0, 0, 0, 4*lambda[0], 4*lambda[1], 4*lambda[2]
418+
};
419+
420+
memcpy(grad, g, sizeof(g));
421+
}
422+
423+
unsigned int cg2_3dshape[] = { 1, 1, 0, 0 };
424+
425+
double cg2_3dnodes[] = { 0, 0, 0,
426+
1, 0, 0,
427+
0, 1, 0,
428+
0, 0, 1,
429+
0.5, 0, 0,
430+
0.5, 0.5, 0,
431+
0, 0.5, 0,
432+
0, 0, 0.5,
433+
0.5, 0, 0.5,
434+
0, 0.5, 0.5 };
435+
436+
eldefninstruction cg2_3deldefn[] = {
437+
LINE(0,0,1), // Identify line subelement with vertex indices (0,1)
438+
LINE(1,1,2), // Identify line subelement with vertex indices (1,2)
439+
LINE(2,2,0), // Identify line subelement with vertex indices (2,0)
440+
LINE(3,0,3), // Identify line subelement with vertex indices (0,3)
441+
LINE(4,1,3), // Identify line subelement with vertex indices (1,3)
442+
LINE(5,2,3), // Identify line subelement with vertex indices (2,3)
443+
QUANTITY(0,0,0), // Fetch quantity on vertex 0
444+
QUANTITY(0,1,0), // Fetch quantity on vertex 1
445+
QUANTITY(0,2,0), // Fetch quantity on vertex 2
446+
QUANTITY(0,3,0), // Fetch quantity on vertex 3
447+
QUANTITY(1,0,0), // Fetch quantity from line 0
448+
QUANTITY(1,1,0), // Fetch quantity from line 1
449+
QUANTITY(1,2,0), // Fetch quantity from line 2
450+
QUANTITY(1,3,0), // Fetch quantity from line 3
451+
QUANTITY(1,4,0), // Fetch quantity from line 4
452+
QUANTITY(1,5,0), // Fetch quantity from line 5
453+
ENDDEFN
454+
};
455+
456+
fespace *cg2_3d_lower[] = {
457+
&cg2_2d,
458+
&cg2_1d,
459+
NULL
460+
};
461+
462+
fespace cg2_3d = {
463+
.name = "CG2",
464+
.grade = 3,
465+
.shape = cg2_3dshape,
466+
.degree = 2,
467+
.nnodes = 10,
468+
.nsubel = 6,
469+
.nodes = cg2_3dnodes,
470+
.ifn = cg2_3dinterpolate,
471+
.gfn = cg2_3dgrad,
472+
.eldefn = cg2_3deldefn,
473+
.lower = cg2_3d_lower
474+
};
475+
476+
/* -------------------------------------------------------
477+
* List of finite elements
478+
* ------------------------------------------------------- */
479+
317480
fespace *fespaces[] = {
318481
&cg1_1d,
319482
&cg2_1d,
320483
&cg1_2d,
321484
&cg2_2d,
485+
&cg1_3d,
486+
&cg2_3d,
322487
NULL
323488
};
324489

‎src/geometry/fespace.h‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,7 @@
77
#ifndef fespace_h
88
#define fespace_h
99

10-
#include "mesh.h"
11-
#include "field.h"
10+
#include "geometry.h"
1211

1312
/* -------------------------------------------------------
1413
* Discretization type definitions

‎src/geometry/functional.c‎

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -4810,15 +4810,7 @@ bool volumeintegral_integrand(vm *v, objectmesh *mesh, elementid id, int nv, int
48104810
/* Set up quantities */
48114811
integral_cleartlvars(v);
48124812
vm_settlvar(v, elementhandle, MORPHO_OBJECT(&elref));
4813-
4814-
value q0[iref.nfields+1], q1[iref.nfields+1], q2[iref.nfields+1], q3[iref.nfields+1];
4815-
value *q[4] = { q0, q1, q2, q3 };
4816-
for (unsigned int k=0; k<iref.nfields; k++) {
4817-
for (unsigned int i=0; i<nv; i++) {
4818-
field_getelement(MORPHO_GETFIELD(iref.fields[k]), MESH_GRADE_VERTEX, vid[i], 0, &q[i][k]);
4819-
}
4820-
}
4821-
4813+
48224814
if (MORPHO_ISDICTIONARY(iref.method)) {
48234815
double err;
48244816
quantity quantities[iref.nfields+1];
@@ -4828,7 +4820,16 @@ bool volumeintegral_integrand(vm *v, objectmesh *mesh, elementid id, int nv, int
48284820
success=integrate(integral_integrandfn, MORPHO_GETDICTIONARY(iref.method), morpho_geterror(v), mesh->dim, MESH_GRADE_VOLUME, x, iref.nfields, quantities, &iref, out, &err);
48294821

48304822
integral_clearquantities(iref.nfields, quantities);
4823+
integral_clearelref(&elref);
48314824
} else {
4825+
value q0[iref.nfields+1], q1[iref.nfields+1], q2[iref.nfields+1], q3[iref.nfields+1];
4826+
value *q[4] = { q0, q1, q2, q3 };
4827+
for (unsigned int k=0; k<iref.nfields; k++) {
4828+
for (unsigned int i=0; i<nv; i++) {
4829+
field_getelement(MORPHO_GETFIELD(iref.fields[k]), MESH_GRADE_VERTEX, vid[i], 0, &q[i][k]);
4830+
}
4831+
}
4832+
48324833
success=integrate_integrate(integral_integrandfn, mesh->dim, MESH_GRADE_VOLUME, x, iref.nfields, q, &iref, out);
48334834
}
48344835

‎src/geometry/integrate.c‎

Lines changed: 84 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1091,6 +1091,10 @@ double tripts[] = {
10911091
0.3333333333333333, 0.3333333333333333, 0.3333333333333333
10921092
};
10931093

1094+
double tri0wts[] = {
1095+
1.0
1096+
};
1097+
10941098
double tri4wts[] = {
10951099
-0.5625, 0.5208333333333332, 0.5208333333333332, 0.5208333333333332
10961100
};
@@ -1141,6 +1145,16 @@ quadraturerule tri4 = {
11411145
.ext = &tri10
11421146
};
11431147

1148+
quadraturerule tri0 = {
1149+
.name = "tri0",
1150+
.grade = 2,
1151+
.order = 1, // 1pt rule is order 1,
1152+
.nnodes = 1,
1153+
.nodes = tripts,
1154+
.weights = tri0wts,
1155+
.ext = &tri4
1156+
};
1157+
11441158
// CUBTRI rule from D. P. Laurie, ACM Transactions on Mathematical Software, Vol 8, No. 2, June 1982,Pages 210-218
11451159

11461160
double cubtripts[] = {
@@ -1166,6 +1180,10 @@ double cubtripts[] = {
11661180
0.232102326775050368,0.738416812340510066,0.0294808608844395667
11671181
};
11681182

1183+
double cubtri0wts[] = {
1184+
1.0
1185+
};
1186+
11691187
double cubtri7wts[] = {
11701188
0.225000000000000000, 0.125939180544827153, 0.125939180544827153,
11711189
0.125939180544827153, 0.132394152788506181, 0.132394152788506181,
@@ -1202,6 +1220,16 @@ quadraturerule cubtri7 = {
12021220
.ext = &cubtri19
12031221
};
12041222

1223+
quadraturerule cubtri0 = {
1224+
.name = "cubtri0",
1225+
.grade = 2,
1226+
.order = 1,
1227+
.nnodes = 1,
1228+
.nodes = cubtripts,
1229+
.weights = cubtri0wts,
1230+
.ext = &cubtri7
1231+
};
1232+
12051233
// Grundmann-Möller embedded rules:
12061234
// SIAM Journal on Numerical Analysis , Apr., 1978, Vol. 15, No. 2 (Apr., 1978), pp. 282-290
12071235
// Computed with simplex_gm_rule [https://people.sc.fsu.edu/~jburkardt/c_src/simplex_gm_rule/simplex_gm_rule.html]
@@ -1273,6 +1301,10 @@ double grundmann2dpts[] = {
12731301
0.07692307692307693, 0.07692307692307693, 0.8461538461538461,
12741302
};
12751303

1304+
double grundmann2d0wts[] = {
1305+
1.0
1306+
};
1307+
12761308
double grundmann2d1wts[] = {
12771309
-0.5625, 0.5208333333333333, 0.5208333333333333, 0.5208333333333333
12781310
};
@@ -1370,6 +1402,16 @@ quadraturerule grundmann2d1 = {
13701402
.ext = &grundmann2d2
13711403
};
13721404

1405+
quadraturerule grundmann2d0 = {
1406+
.name = "grundmann2d0",
1407+
.grade = 2,
1408+
.order = 1,
1409+
.nnodes = 1,
1410+
.nodes = grundmann2dpts,
1411+
.weights = grundmann2d0wts,
1412+
.ext = &grundmann2d1
1413+
};
1414+
13731415
/* --------------------------------
13741416
* Tetrahedron
13751417
* -------------------------------- */
@@ -1746,6 +1788,10 @@ double grundmann3dpts[] = {
17461788
0.071428571428571428571,0.071428571428571428571,0.071428571428571428571
17471789
};
17481790

1791+
double grundmann3d0wts[] = {
1792+
1.0
1793+
};
1794+
17491795
double grundmann3d1wts[] = {
17501796
-0.8,0.45,0.45,0.45,0.45
17511797
};
@@ -1895,6 +1941,16 @@ quadraturerule grundmann3d1 = {
18951941
.ext = &grundmann3d2
18961942
};
18971943

1944+
quadraturerule grundmann3d0 = {
1945+
.name = "grundmann3d0",
1946+
.grade = 3,
1947+
.order = 1,
1948+
.nnodes = 1,
1949+
.nodes = grundmann3dpts,
1950+
.weights = grundmann3d0wts,
1951+
.ext = &grundmann3d1
1952+
};
1953+
18981954
/* --------------------------------
18991955
* List of quadrature rules
19001956
* -------------------------------- */
@@ -1906,14 +1962,22 @@ quadraturerule *quadrules[] = {
19061962
&gauss5, &kronrod11,
19071963
&gauss7, &kronrod15,
19081964

1909-
&tri4, &tri10, &tri20,
1910-
&cubtri7, &cubtri19,
1911-
&grundmann2d1, &grundmann2d2, &grundmann2d3, &grundmann2d4, &grundmann2d5,
1965+
&tri0, &tri4, &tri10, &tri20,
1966+
&cubtri0, &cubtri7, &cubtri19,
1967+
&grundmann2d0, &grundmann2d1, &grundmann2d2, &grundmann2d3, &grundmann2d4, &grundmann2d5,
19121968

19131969
&keast4, &keast5,
19141970
&tet5, &tet6,
19151971

1916-
&grundmann3d1, &grundmann3d2, &grundmann3d3, &grundmann3d4, &grundmann3d5,
1972+
&grundmann3d0, &grundmann3d1, &grundmann3d2, &grundmann3d3, &grundmann3d4, &grundmann3d5,
1973+
NULL
1974+
};
1975+
1976+
// Specify a list of default rules for each grade
1977+
quadraturerule *defaultquadrule[] = {
1978+
&gauss5,
1979+
&cubtri7,
1980+
&grundmann3d0,
19171981
NULL
19181982
};
19191983

@@ -2411,7 +2475,8 @@ bool integrator_quadrature(integrator *integrate, quadraturerule *rule, quadratu
24112475
eps[ip]=fabs(r[ip]-r[ip-1]);
24122476
nmin = q->nnodes;
24132477

2414-
if (fabs(eps[ip]/r[ip])<1e-14) break;
2478+
if (fabs(r[ip])<integrate->ztol ||
2479+
fabs(eps[ip]/r[ip])<integrate->tol) break;
24152480
}
24162481

24172482
work->lval = work->weight*r[ip-1];
@@ -2560,6 +2625,17 @@ bool integrator_matchrulebyorder(int grade, int minorder, int maxorder, bool hig
25602625
return (best>=0);
25612626
}
25622627

2628+
/** Returns a default rule for each grade */
2629+
bool integrator_matchrulebygrade(int grade, quadraturerule **out) {
2630+
for (int i=0; defaultquadrule[i]!=NULL; i++) {
2631+
if (defaultquadrule[i]->grade==grade) {
2632+
*out = defaultquadrule[i];
2633+
return true;
2634+
}
2635+
}
2636+
return false;
2637+
}
2638+
25632639
/** Configures an integrator based on the grade to integrate and hints for order and rule type
25642640
* @param[in] integrate - integrator structure to be configured
25652641
* @param[in] err - error structure to report errors to
@@ -2580,8 +2656,10 @@ bool integrator_configure(integrator *integrate, error *err, bool adapt, int gra
25802656
error_writewithid(err, INTEGRATE_RLNTFND, name);
25812657
return false;
25822658
}
2659+
} else if (order>=0) {
2660+
integrator_matchrulebyorder(grade, order, INT_MAX, false, &integrate->rule);
25832661
} else {
2584-
integrator_matchrulebyorder(grade, (order<0 ? 0 : order), INT_MAX, (order<0), &integrate->rule);
2662+
integrator_matchrulebygrade(grade, &integrate->rule);
25852663
}
25862664

25872665
// Check we succeeded in finding a rule

0 commit comments

Comments
 (0)