1
#ifdef PETSC_RCS_HEADER
2
static char vcid[] = "$Id: quadratic.c,v 1.7 2000/01/10 03:54:16 knepley Exp $";
6
Defines piecewise quadratic function space on a two dimensional
7
grid. Suitable for finite element type discretization of a PDE.
10
#include "src/grid/discretization/discimpl.h" /*I "discretization.h" I*/
11
#include "src/mesh/impls/triangular/triimpl.h"
13
/* For precomputed integrals, the table is structured as follows:
15
precompInt[op,i,j] = \int_{SE} <op \phi^i(\xi,\eta), \phi^j(\xi,\eta)> |J^{-1}|
17
The Jacobian in this case may not be constant over the element in question.
18
The numbering of the nodes in the standard element is
30
#define __FUNCT__ "DiscDestroy_Triangular_2D_Quadratic"
31
static int DiscDestroy_Triangular_2D_Quadratic(Discretization disc) {
33
PetscFunctionReturn(0);
37
#define __FUNCT__ "DiscView_Triangular_2D_Quadratic_File"
38
static int DiscView_Triangular_2D_Quadratic_File(Discretization disc, PetscViewer viewer) {
40
PetscViewerASCIIPrintf(viewer, "Quadratic discretization\n");
41
PetscViewerASCIIPrintf(viewer, " %d shape functions per component\n", disc->funcs);
42
PetscViewerASCIIPrintf(viewer, " %d registered operators\n", disc->numOps);
43
PetscFunctionReturn(0);
47
#define __FUNCT__ "DiscView_Triangular_2D_Quadratic"
48
static int DiscView_Triangular_2D_Quadratic(Discretization disc, PetscViewer viewer) {
53
ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &isascii); CHKERRQ(ierr);
54
if (isascii == PETSC_TRUE) {
55
ierr = DiscView_Triangular_2D_Quadratic_File(disc, viewer); CHKERRQ(ierr);
57
PetscFunctionReturn(0);
61
#define __FUNCT__ "DiscEvaluateShapeFunctions_Triangular_2D_Quadratic_Private"
62
int DiscEvaluateShapeFunctions_Triangular_2D_Quadratic_Private(double xi, double eta, double *coords, double *x,
63
double *y, double *dxxi, double *dxet, double *dyxi, double *dyet)
65
/* ASSUMPTION: The coordinates passed in are corrected for periodicity */
68
*dxxi = 0.0; *dxet = 0.0;
69
*dyxi = 0.0; *dyet = 0.0;
70
/* \phi^0: 1 - 3 (\xi + \eta) + 2 (\xi + \eta)^2 */
71
*x += coords[0*2+0]*(1.0 - (xi + eta))*(1.0 - 2.0*(xi + eta));
72
*dxxi += coords[0*2+0]*(-3.0 + 4.0*(xi + eta));
73
*dxet += coords[0*2+0]*(-3.0 + 4.0*(xi + eta));
74
*y += coords[0*2+1]*(1.0 - (xi + eta))*(1.0 - 2.0*(xi + eta));
75
*dyxi += coords[0*2+1]*(-3.0 + 4.0*(xi + eta));
76
*dyet += coords[0*2+1]*(-3.0 + 4.0*(xi + eta));
77
/* \phi^1: \xi (2\xi - 1) */
78
*x += coords[1*2+0]*(xi*(2.0*xi - 1.0));
79
*dxxi += coords[1*2+0]*(4.0*xi - 1.0);
81
*y += coords[1*2+1]*(xi*(2.0*xi - 1.0));
82
*dyxi += coords[1*2+1]*(4.0*xi - 1.0);
84
/* \phi^2: \eta (2\eta - 1) */
85
*x += coords[2*2+0]*(eta*(2.0*eta - 1.0));
87
*dxet += coords[2*2+0]*(4.0*eta - 1.0);
88
*y += coords[2*2+1]*(eta*(2.0*eta - 1.0));
90
*dyet += coords[2*2+1]*(4.0*eta - 1.0);
91
/* \phi^3: 4 \xi \eta */
92
*x += coords[3*2+0]*(4.0*xi*eta);
93
*dxxi += coords[3*2+0]*(4.0*eta);
94
*dxet += coords[3*2+0]*(4.0*xi);
95
*y += coords[3*2+1]*(4.0*xi*eta);
96
*dyxi += coords[3*2+1]*(4.0*eta);
97
*dyet += coords[3*2+1]*(4.0*xi);
98
/* \phi^4: 4 \eta (1 - \xi - \eta) */
99
*x += coords[4*2+0]*(4.0*eta*(1.0 - (xi + eta)));
100
*dxxi += coords[4*2+0]*(-4.0*eta);
101
*dxet += coords[4*2+0]*(-8.0*eta + 4.0*(1.0 - xi));
102
*y += coords[4*2+1]*(4.0*eta*(1.0 - (xi + eta)));
103
*dyxi += coords[4*2+1]*(-4.0*eta);
104
*dyet += coords[4*2+1]*(-8.0*eta + 4.0*(1.0 - xi));
105
/* \phi^5: 4 \xi (1 - \xi - \eta) */
106
*x += coords[5*2+0]*(4.0*xi*(1.0 - (xi + eta)));
107
*dxxi += coords[5*2+0]*(-8.0*xi + 4.0*(1.0 - eta));
108
*dxet += coords[5*2+0]*(-4.0*xi);
109
*y += coords[5*2+1]*(4.0*xi*(1.0 - (xi + eta)));
110
*dyxi += coords[5*2+1]*(-8.0*xi + 4.0*(1.0 - eta));
111
*dyet += coords[5*2+1]*(-4.0*xi);
112
PetscLogFlops(36+20+20+20+30+20);
113
PetscFunctionReturn(0);
117
#define __FUNCT__ "DiscTransformCoords_Triangular_2D_Quadratic"
118
int DiscTransformCoords_Triangular_2D_Quadratic(double x, double y, double *coords, double *newXi, double *newEta)
120
/* ASSUMPTION: The coordinates passed in are corrected for periodicity */
121
double xi, eta; /* Canonical coordinates of the point */
122
double dxix; /* \PartDer{\xi}{x} */
123
double detx; /* \PartDer{\eta}{x} */
124
double dxiy; /* \PartDer{\xi}{y} */
125
double dety; /* \PartDer{\eta}{y} */
126
double dxxi; /* \PartDer{x}{\xi} */
127
double dxet; /* \PartDer{x}{\eta} */
128
double dyxi; /* \PartDer{y}{\xi} */
129
double dyet; /* \PartDer{y}{\eta} */
130
double new_x; /* x(\xi,\eta) */
131
double new_y; /* x(\xi,\eta) */
132
double residual_x; /* x(\xi,\eta) - x */
133
double residual_y; /* x(\xi,\eta) - y */
134
double jac, invjac; /* The Jacobian determinant and its inverse */
142
\sum_f x(f) \phi^f(\xi,\eta) = x
143
\sum_f y(f) \phi^f(\xi,\eta) = y
145
where f runs over nodes (each one has coordinates and a shape function). We
146
will use Newton's method
148
/ \sum_f x(f) \PartDer{\phi^f}{\xi} \sum_f x(f) \PartDer{\phi^f}{\eta} \ / d\xi \ = / \sum_f x(f) \phi^f(\xi,\eta) - x\
149
\ \sum_f y(f) \PartDer{\phi^f}{\xi} \sum_f y(f) \PartDer{\phi^f}{\eta} / \ d\eta / \ \sum_f y(f) \phi^f(\xi,\eta) - y/
151
which can be rewritten more compactly as
153
/ \PartDer{x}{\xi} \PartDer{x}{\eta} \ / d\xi \ = / x(\xi,\eta) - x \
154
\ \PartDer{y}{\xi} \PartDer{y}{\eta} / \ d\eta / \ y(\xi,\eta) - y /
156
The initial guess will be the linear solution.
158
ASSUMPTION: The coordinates passed in are all on the same sheet as x,y
160
/* Form linear solution */
161
dxxi = coords[1*2+0] - coords[0*2+0];
162
dxet = coords[2*2+0] - coords[0*2+0];
163
dyxi = coords[1*2+1] - coords[0*2+1];
164
dyet = coords[2*2+1] - coords[0*2+1];
165
jac = PetscAbsReal(dxxi*dyet - dxet*dyxi);
166
if (jac < 1.0e-14) SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
172
xi = dxix*(x - coords[0*2+0]) + dxiy*(y - coords[0*2+1]);
173
eta = detx*(x - coords[0*2+0]) + dety*(y - coords[0*2+1]);
174
for(iter = 0; iter < maxIters; iter++) {
175
/* This is clumsy, but I can't think of anything better right now */
176
ierr = DiscEvaluateShapeFunctions_Triangular_2D_Quadratic_Private(xi, eta, coords, &new_x, &new_y, &dxxi, &dxet, &dyxi, &dyet);
179
/* Check for convergence -- I should maybe make the tolerance variable */
180
residual_x = new_x - x;
181
residual_y = new_y - y;
182
if (PetscAbsReal(residual_x) + PetscAbsReal(residual_y) < 1.0e-6) break;
184
/* Solve the system */
185
jac = PetscAbsReal(dxxi*dyet - dxet*dyxi);
191
/* These are the elements of the inverse matrix */
197
xi -= dxix*residual_x + dxiy*residual_y;
198
eta -= detx*residual_x + dety*residual_y;
200
if (iter == maxIters) {
201
PetscLogInfo(PETSC_NULL, "DiscTransformCoords_Triangular_2D_Quadratic: Newton iteration did not converge\n");
202
PetscLogInfo(PETSC_NULL, "x: %g y: %g maxIters: %d\n", x, y, maxIters);
203
for(iter = 0; iter < 6; iter++) {
204
PetscLogInfo(PETSC_NULL, " x%d: %g y%d: %g\n", iter, coords[iter*2+0], iter, coords[iter*2+1]);
206
/* Use linear interpolation */
207
xi = dxix*(x - coords[0*2+0]) + dxiy*(y - coords[0*2+1]);
208
eta = detx*(x - coords[0*2+0]) + dety*(y - coords[0*2+1]);
213
PetscLogFlops(7+15+19*iter);
214
PetscFunctionReturn(0);
218
#define __FUNCT__ "DiscEvaluateFunctionGalerkin_Triangular_2D_Quadratic"
219
static int DiscEvaluateFunctionGalerkin_Triangular_2D_Quadratic(Discretization disc, Mesh mesh, PointFunction f,
220
PetscScalar alpha, int elem, PetscScalar *array, void *ctx)
222
Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
223
double *nodes = tri->nodes;
224
int *elements = tri->faces;
225
int numCorners = mesh->numCorners;
226
int comp = disc->comp; /* The number of components in this field */
227
int funcs = disc->funcs; /* The number of shape functions per component */
228
PetscScalar *funcVal = disc->funcVal; /* Function value at a quadrature point */
229
int numQuadPoints = disc->numQuadPoints; /* Number of points used for Gaussian quadrature */
230
double *quadWeights = disc->quadWeights; /* Weights in the standard element for Gaussian quadrature */
231
double *quadShapeFuncs = disc->quadShapeFuncs; /* Shape functions evaluated at quadrature points */
232
double *quadShapeFuncDers = disc->quadShapeFuncDers; /* Shape function derivatives at quadrature points */
233
double jac; /* |J| for map to standard element */
234
double x, y; /* The integration point */
235
double dxxi; /* \PartDer{x}{\xi} */
236
double dxet; /* \PartDer{x}{\xi} */
237
double dyxi; /* \PartDer{y}{\eta} */
238
double dyet; /* \PartDer{y}{\eta} */
239
double coords[MAX_CORNERS*2];
241
int i, j, k, func, p;
242
#ifdef PETSC_USE_BOPT_g
248
ierr = MPI_Comm_rank(disc->comm, &rank); CHKERRQ(ierr);
250
/* For dummy collective calls */
251
if (array == PETSC_NULL) {
252
for(p = 0; p < numQuadPoints; p++) {
253
ierr = (*f)(0, 0, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, ctx); CHKERRQ(ierr);
255
PetscFunctionReturn(0);
258
#ifdef PETSC_USE_BOPT_g
259
if ((elem < 0) || (elem >= mesh->part->numOverlapElements)) {
260
SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid element %d should be in [0,%d)", elem, mesh->part->numOverlapElements);
263
/* Calculate the determinant of the inverse Jacobian of the map to the standard element */
264
for(i = 0; i < numCorners; i++) {
265
coords[i*2] = nodes[elements[elem*numCorners+i]*2];
266
coords[i*2+1] = nodes[elements[elem*numCorners+i]*2+1];
269
/* Check for constant jacobian here */
271
jac = PetscAbsReal((coords[2] - coords[0])*(coords[5] - coords[1]) - (coords[4] - coords[0])*(coords[3] - coords[1]));
273
PetscPrintf(PETSC_COMM_SELF, "[%d], elem: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
274
rank, elem, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
275
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
278
#ifdef PETSC_USE_BOPT_g
279
ierr = PetscOptionsHasName(PETSC_NULL, "-trace_assembly", &opt); CHKERRQ(ierr);
280
if (opt == PETSC_TRUE) {
281
PetscPrintf(PETSC_COMM_SELF, "[%d]elem: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
282
rank, elem, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
283
PetscPrintf(PETSC_COMM_SELF, " x4: %g y4: %g x5: %g y5: %g x6: %g y6: %g\n",
284
coords[6], coords[7], coords[8], coords[9], coords[10], coords[11]);
288
/* Calculate element vector entries by Gaussian quadrature */
289
for(p = 0; p < numQuadPoints; p++)
291
/* x = \sum^{funcs}_{f=1} x_f \phi^f(p)
292
\PartDer{x}{\xi}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\xi}
293
\PartDer{x}{\eta}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\eta}
294
y = \sum^{funcs}_{f=1} y_f \phi^f(p)
295
\PartDer{y}{\xi}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\xi}
296
\PartDer{y}{\eta}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\eta} */
298
dxxi = 0.0; dxet = 0.0;
299
dyxi = 0.0; dyet = 0.0;
300
if (mesh->isPeriodic == PETSC_TRUE) {
301
for(func = 0; func < funcs; func++)
303
x += MeshPeriodicRelativeX(mesh, coords[func*2], coords[0])*quadShapeFuncs[p*funcs+func];
304
dxxi += MeshPeriodicRelativeX(mesh, coords[func*2], coords[0])*quadShapeFuncDers[p*funcs*2+func*2];
305
dxet += MeshPeriodicRelativeX(mesh, coords[func*2], coords[0])*quadShapeFuncDers[p*funcs*2+func*2+1];
306
y += MeshPeriodicRelativeY(mesh, coords[func*2+1], coords[1])*quadShapeFuncs[p*funcs+func];
307
dyxi += MeshPeriodicRelativeY(mesh, coords[func*2+1], coords[1])*quadShapeFuncDers[p*funcs*2+func*2];
308
dyet += MeshPeriodicRelativeY(mesh, coords[func*2+1], coords[1])*quadShapeFuncDers[p*funcs*2+func*2+1];
310
x = MeshPeriodicX(mesh, x);
311
y = MeshPeriodicY(mesh, y);
313
for(func = 0; func < funcs; func++)
315
x += coords[func*2] *quadShapeFuncs[p*funcs+func];
316
dxxi += coords[func*2] *quadShapeFuncDers[p*funcs*2+func*2];
317
dxet += coords[func*2] *quadShapeFuncDers[p*funcs*2+func*2+1];
318
y += coords[func*2+1]*quadShapeFuncs[p*funcs+func];
319
dyxi += coords[func*2+1]*quadShapeFuncDers[p*funcs*2+func*2];
320
dyet += coords[func*2+1]*quadShapeFuncDers[p*funcs*2+func*2+1];
323
jac = PetscAbsReal(dxxi*dyet - dxet*dyxi);
325
PetscPrintf(PETSC_COMM_SELF, "[%d]p: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
326
rank, p, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
327
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
329
ierr = (*f)(1, comp, &x, &y, PETSC_NULL, funcVal, ctx); CHKERRQ(ierr);
330
#ifdef PETSC_USE_BOPT_g
331
ierr = PetscOptionsHasName(PETSC_NULL, "-trace_assembly", &opt); CHKERRQ(ierr);
332
if (opt == PETSC_TRUE) {
333
PetscPrintf(PETSC_COMM_SELF, "[%d]p: %d jac: %g", rank, p, jac);
334
for(j = 0; j < comp; j++)
335
PetscPrintf(PETSC_COMM_SELF, " func[%d]: %g", j, PetscRealPart(funcVal[j]));
336
PetscPrintf(PETSC_COMM_SELF, "\n");
340
for(i = 0, k = 0; i < funcs; i++) {
341
for(j = 0; j < comp; j++, k++) {
342
array[k] += alpha*funcVal[j]*quadShapeFuncs[p*funcs+i]*jac*quadWeights[p];
343
#ifdef PETSC_USE_BOPT_g
344
ierr = PetscOptionsHasName(PETSC_NULL, "-trace_assembly", &opt); CHKERRQ(ierr);
345
if (opt == PETSC_TRUE) {
346
PetscPrintf(PETSC_COMM_SELF, "[%d] array[%d]: %g\n", rank, k, PetscRealPart(array[k]));
352
PetscLogFlops((3 + 12*funcs + 5*funcs*comp) * numQuadPoints);
353
PetscFunctionReturn(0);
357
#define __FUNCT__ "DiscEvaluateOperatorGalerkin_Triangular_2D_Quadratic"
358
static int DiscEvaluateOperatorGalerkin_Triangular_2D_Quadratic(Discretization disc, Mesh mesh, int elemSize,
359
int rowStart, int colStart, int op, PetscScalar alpha,
360
int elem, PetscScalar *field, PetscScalar *array, void *ctx)
362
Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
363
double *nodes = tri->nodes; /* The node coordinates */
364
int *elements = tri->faces; /* The element corners */
365
int numCorners = mesh->numCorners; /* The number of corners per element */
366
Operator oper = disc->operators[op]; /* The operator to discretize */
367
Discretization test = oper->test; /* The space of test functions */
368
OperatorFunction opFunc = oper->opFunc; /* Integrals of operators which depend on J */
369
PetscScalar *precompInt = oper->precompInt; /* Precomputed integrals of the operator on shape functions */
370
int rowSize = test->size; /* The number of shape functions per element */
371
int colSize = disc->size; /* The number of test functions per element */
372
double x21, x31, y21, y31; /* Coordinates of the element, with point 1 at the origin */
373
double jac; /* |J| for map to standard element */
374
double coords[MAX_CORNERS*2]; /* Coordinates of the element */
380
ierr = MPI_Comm_rank(disc->comm, &rank); CHKERRQ(ierr);
381
#ifdef PETSC_USE_BOPT_g
382
/* Check for valid operator */
383
if ((op < 0) || (op >= disc->numOps) || (!disc->operators[op])) SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid operator");
386
if (precompInt != PETSC_NULL)
388
/* Calculate the determinant of the inverse Jacobian of the map to the standard element
389
which has been specified as constant here - 1/|x_{21} y_{31} - x_{31} y_{21}| */
390
if (mesh->isPeriodic == PETSC_TRUE) {
391
x21 = MeshPeriodicDiffX(mesh, nodes[elements[elem*numCorners+1]*2] - nodes[elements[elem*numCorners]*2]);
392
x31 = MeshPeriodicDiffX(mesh, nodes[elements[elem*numCorners+2]*2] - nodes[elements[elem*numCorners]*2]);
393
y21 = MeshPeriodicDiffY(mesh, nodes[elements[elem*numCorners+1]*2+1] - nodes[elements[elem*numCorners]*2+1]);
394
y31 = MeshPeriodicDiffY(mesh, nodes[elements[elem*numCorners+2]*2+1] - nodes[elements[elem*numCorners]*2+1]);
396
x21 = nodes[elements[elem*numCorners+1]*2] - nodes[elements[elem*numCorners]*2];
397
x31 = nodes[elements[elem*numCorners+2]*2] - nodes[elements[elem*numCorners]*2];
398
y21 = nodes[elements[elem*numCorners+1]*2+1] - nodes[elements[elem*numCorners]*2+1];
399
y31 = nodes[elements[elem*numCorners+2]*2+1] - nodes[elements[elem*numCorners]*2+1];
401
jac = PetscAbsReal(x21*y31 - x31*y21);
403
PetscPrintf(PETSC_COMM_SELF, "[%d]x21: %g y21: %g x31: %g y31: %g jac: %g\n", rank, x21, y21, x31, y31, jac);
404
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
406
/* PetscPrintf(PETSC_COMM_SELF, "x21: %g y21: %g x31: %g y31: %g\n", x21, y21, x31, y31, jac); */
408
/* Calculate element matrix entries which are all precomputed */
409
for(i = 0; i < rowSize; i++)
410
for(j = 0; j < colSize; j++)
411
array[(i+rowStart)*elemSize + j+colStart] += alpha*precompInt[i*colSize + j]*jac;
412
PetscLogFlops(7 + 2*rowSize*colSize);
416
if (opFunc == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid function");
417
if (mesh->isPeriodic == PETSC_TRUE) {
418
coords[0*2+0] = nodes[elements[elem*numCorners+0]*2+0];
419
coords[0*2+1] = nodes[elements[elem*numCorners+0]*2+1];
420
for(f = 1; f < PetscMax(disc->funcs, test->funcs); f++) {
421
coords[f*2+0] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+f]*2+0], coords[0*2+0]);
422
coords[f*2+1] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+f]*2+1], coords[0*2+1]);
425
for(f = 0; f < PetscMax(disc->funcs, test->funcs); f++) {
426
coords[f*2+0] = nodes[elements[elem*numCorners+f]*2+0];
427
coords[f*2+1] = nodes[elements[elem*numCorners+f]*2+1];
431
ierr = (*opFunc)(disc, test, rowSize, colSize, rowStart, colStart, elemSize, coords, alpha, field, array, ctx);
434
PetscFunctionReturn(0);
438
#define __FUNCT__ "DiscEvaluateNonlinearOperatorGalerkin_Triangular_2D_Quadratic"
439
static int DiscEvaluateNonlinearOperatorGalerkin_Triangular_2D_Quadratic(Discretization disc, Mesh mesh, NonlinearOperator f,
440
PetscScalar alpha, int elem, int numArgs, PetscScalar **field,
441
PetscScalar *vec, void *ctx)
443
Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
444
double *nodes = tri->nodes;
445
int *elements = tri->faces;
446
int numCorners = mesh->numCorners;
447
int comp = disc->comp; /* The number of components in this field */
448
int funcs = disc->funcs; /* The number of shape functions per component */
449
PetscScalar *funcVal = disc->funcVal; /* Function value at a quadrature point */
450
PetscScalar **fieldVal = disc->fieldVal; /* Field value and derivatives at a quadrature point */
451
double jac; /* |J| for map to standard element */
452
double invjac; /* |J^{-1}| for map from standard element */
453
int numQuadPoints; /* Number of points used for Gaussian quadrature */
454
double *quadWeights; /* Weights in the standard element for Gaussian quadrature */
455
double *quadShapeFuncs; /* Shape functions evaluated at quadrature points */
456
double *quadShapeFuncDers; /* Shape function derivatives evaluated at quadrature points */
457
double x, y; /* The integration point */
458
double dxxi; /* \PartDer{x}{\xi} */
459
double dxet; /* \PartDer{x}{\eta} */
460
double dyxi; /* \PartDer{y}{\xi} */
461
double dyet; /* \PartDer{y}{\eta} */
462
double dxix; /* \PartDer{\xi}{x} */
463
double detx; /* \PartDer{\eta}{x} */
464
double dxiy; /* \PartDer{\xi}{y} */
465
double dety; /* \PartDer{\eta}{y} */
466
PetscScalar dfxi; /* \PartDer{field}{\xi} */
467
PetscScalar dfet; /* \PartDer{field}{\eta} */
468
double coords[12]; /* Coordinates of the element */
470
int i, j, k, func, p, arg;
471
#ifdef PETSC_USE_BOPT_g
477
ierr = MPI_Comm_rank(disc->comm, &rank); CHKERRQ(ierr);
478
numQuadPoints = disc->numQuadPoints;
479
quadWeights = disc->quadWeights;
480
quadShapeFuncs = disc->quadShapeFuncs;
481
quadShapeFuncDers = disc->quadShapeFuncDers;
483
/* Calculate the determinant of the inverse Jacobian of the map to the standard element */
484
if (mesh->isPeriodic == PETSC_TRUE) {
485
coords[0*2+0] = nodes[elements[elem*numCorners+0]*2+0];
486
coords[0*2+1] = nodes[elements[elem*numCorners+0]*2+1];
487
for(func = 1; func < funcs; func++) {
488
coords[func*2+0] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+func]*2+0], coords[0*2+0]);
489
coords[func*2+1] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+func]*2+1], coords[0*2+1]);
492
for(func = 0; func < funcs; func++) {
493
coords[func*2+0] = nodes[elements[elem*numCorners+func]*2+0];
494
coords[func*2+1] = nodes[elements[elem*numCorners+func]*2+1];
497
/* Check for constant jacobian here */
499
jac = PetscAbsReal((coords[2] - coords[0])*(coords[5] - coords[1]) - (coords[4] - coords[0])*(coords[3] - coords[1]));
501
PetscPrintf(PETSC_COMM_SELF, "[%d], elem: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
502
rank, elem, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
503
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
506
#ifdef PETSC_USE_BOPT_g
507
ierr = PetscOptionsHasName(PETSC_NULL, "-trace_assembly", &opt); CHKERRQ(ierr);
508
if (opt == PETSC_TRUE) {
509
PetscPrintf(PETSC_COMM_SELF, "[%d]elem: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
510
rank, elem, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
511
PetscPrintf(PETSC_COMM_SELF, " x4: %g y4: %g x5: %g y5: %g x6: %g y6: %g\n",
512
coords[6], coords[7], coords[8], coords[9], coords[10], coords[11]);
516
/* Calculate element vector entries by Gaussian quadrature */
517
for(p = 0; p < numQuadPoints; p++) {
518
/* x = \sum^{funcs}_{f=1} x_f \phi^f(p)
519
\PartDer{x}{\xi}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\xi}
520
\PartDer{x}{\eta}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\eta}
521
y = \sum^{funcs}_{f=1} y_f \phi^f(p)
522
\PartDer{y}{\xi}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\xi}
523
\PartDer{y}{\eta}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\eta}
524
u^i = \sum^{funcs}_{f=1} u^i_f \phi^f(p)
525
\PartDer{u^i}{\xi}(p) = \sum^{funcs}_{f=1} u^i_f \PartDer{\phi^f(p)}{\xi}
526
\PartDer{u^i}{\eta}(p) = \sum^{funcs}_{f=1} u^i_f \PartDer{\phi^f(p)}{\eta} */
528
dxxi = 0.0; dyxi = 0.0;
529
dxet = 0.0; dyet = 0.0;
530
for(arg = 0; arg < numArgs; arg++)
531
for(j = 0; j < comp*3; j++)
532
fieldVal[arg][j] = 0.0;
533
for(func = 0; func < funcs; func++) {
534
x += coords[func*2] *quadShapeFuncs[p*funcs+func];
535
dxxi += coords[func*2] *quadShapeFuncDers[p*funcs*2+func*2];
536
dxet += coords[func*2] *quadShapeFuncDers[p*funcs*2+func*2+1];
537
y += coords[func*2+1]*quadShapeFuncs[p*funcs+func];
538
dyxi += coords[func*2+1]*quadShapeFuncDers[p*funcs*2+func*2];
539
dyet += coords[func*2+1]*quadShapeFuncDers[p*funcs*2+func*2+1];
540
for(arg = 0; arg < numArgs; arg++) {
541
for(j = 0; j < comp; j++) {
542
fieldVal[arg][j*3] += field[arg][func*comp+j]*quadShapeFuncs[p*funcs+func];
543
fieldVal[arg][j*3+1] += field[arg][func*comp+j]*quadShapeFuncDers[p*funcs*2+func*2];
544
fieldVal[arg][j*3+2] += field[arg][func*comp+j]*quadShapeFuncDers[p*funcs*2+func*2+1];
548
if (mesh->isPeriodic == PETSC_TRUE) {
549
x = MeshPeriodicX(mesh, x);
550
y = MeshPeriodicY(mesh, y);
552
jac = PetscAbsReal(dxxi*dyet - dxet*dyxi);
554
PetscPrintf(PETSC_COMM_SELF, "[%d]p: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
555
rank, p, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
556
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
558
/* These are the elements of the inverse matrix */
565
/* Convert the field derivatives to old coordinates */
566
for(arg = 0; arg < numArgs; arg++)
567
for(j = 0; j < comp; j++) {
568
dfxi = fieldVal[arg][j*3+1];
569
dfet = fieldVal[arg][j*3+2];
570
fieldVal[arg][j*3+1] = dfxi*dxix + dfet*detx;
571
fieldVal[arg][j*3+2] = dfxi*dxiy + dfet*dety;
574
ierr = (*f)(1, comp, &x, &y, PETSC_NULL, numArgs, fieldVal, funcVal, ctx); CHKERRQ(ierr);
575
#ifdef PETSC_USE_BOPT_g
576
ierr = PetscOptionsHasName(PETSC_NULL, "-trace_assembly", &opt); CHKERRQ(ierr);
577
if (opt == PETSC_TRUE) {
578
PetscPrintf(PETSC_COMM_SELF, "[%d]p: %d jac: %g", rank, p, jac);
579
for(j = 0; j < comp; j++)
580
PetscPrintf(PETSC_COMM_SELF, " func[%d]: %g", j, PetscRealPart(funcVal[j]));
581
PetscPrintf(PETSC_COMM_SELF, "\n");
585
for(i = 0, k = 0; i < funcs; i++) {
586
for(j = 0; j < comp; j++, k++) {
587
vec[k] += alpha*funcVal[j]*quadShapeFuncs[p*funcs+i]*jac*quadWeights[p];
588
#ifdef PETSC_USE_BOPT_g
589
ierr = PetscOptionsHasName(PETSC_NULL, "-trace_assembly", &opt); CHKERRQ(ierr);
590
if (opt == PETSC_TRUE) {
591
PetscPrintf(PETSC_COMM_SELF, "[%d] vec[%d]: %g\n", rank, k, PetscRealPart(vec[k]));
597
PetscLogFlops(((12 + (6*numArgs + 5)*comp)*funcs + 8 + 6*numArgs*comp) * numQuadPoints);
598
PetscFunctionReturn(0);
602
#define __FUNCT__ "DiscEvaluateALEOperatorGalerkin_Triangular_2D_Quadratic"
603
static int DiscEvaluateALEOperatorGalerkin_Triangular_2D_Quadratic(Discretization disc, Mesh mesh, int elemSize,
604
int rowStart, int colStart, int op, PetscScalar alpha,
605
int elem, PetscScalar *field, PetscScalar *ALEfield, PetscScalar *array,
608
Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
609
double *nodes = tri->nodes;
610
int *elements = tri->faces;
611
int numCorners = mesh->numCorners;
612
Discretization test; /* The space of test functions */
613
Operator oper; /* The operator to discretize */
614
int rowSize; /* The number of shape functions per element */
615
int colSize; /* The number of test functions per element */
616
ALEOperatorFunction opFunc; /* Integrals of operators which depend on J */
617
double coords[MAX_CORNERS*2]; /* Coordinates of the element */
622
#ifdef PETSC_USE_BOPT_g
623
/* Check for valid operator */
624
if ((op < 0) || (op >= disc->numOps) || (!disc->operators[op])) SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid operator");
626
/* Get discretization info */
627
oper = disc->operators[op];
628
opFunc = oper->ALEOpFunc;
630
rowSize = test->size;
631
colSize = disc->size;
633
if (opFunc == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid function");
634
if (mesh->isPeriodic == PETSC_TRUE) {
635
coords[0*2+0] = nodes[elements[elem*numCorners+0]*2+0];
636
coords[0*2+1] = nodes[elements[elem*numCorners+0]*2+1];
637
for(f = 1; f < PetscMax(disc->funcs, test->funcs); f++) {
638
coords[f*2+0] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+f]*2+0], coords[0*2+0]);
639
coords[f*2+1] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+f]*2+1], coords[0*2+1]);
642
for(f = 0; f < PetscMax(disc->funcs, test->funcs); f++) {
643
coords[f*2+0] = nodes[elements[elem*numCorners+f]*2+0];
644
coords[f*2+1] = nodes[elements[elem*numCorners+f]*2+1];
648
ierr = (*opFunc)(disc, test, rowSize, colSize, rowStart, colStart, elemSize, coords, alpha, field, ALEfield, array, ctx);
650
PetscFunctionReturn(0);
654
#define __FUNCT__ "DiscEvaluateNonlinearALEOperatorGalerkin_Triangular_2D_Quadratic"
655
static int DiscEvaluateNonlinearALEOperatorGalerkin_Triangular_2D_Quadratic(Discretization disc, Mesh mesh, NonlinearOperator f,
656
PetscScalar alpha, int elem, int numArgs, PetscScalar **field,
657
PetscScalar *ALEfield, PetscScalar *vec, void *ctx)
659
Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
660
double *nodes = tri->nodes;
661
int *elements = tri->faces;
662
int numCorners = mesh->numCorners;
663
int comp = disc->comp; /* The number of components in this field */
664
int funcs = disc->funcs; /* The number of shape functions per component */
665
PetscScalar *funcVal = disc->funcVal; /* Function value at a quadrature point */
666
PetscScalar **fieldVal = disc->fieldVal; /* Field value and derivatives at a quadrature point */
667
double jac; /* |J| for map to standard element */
668
double invjac; /* |J^{-1}| for map from standard element */
669
int numQuadPoints; /* Number of points used for Gaussian quadrature */
670
double *quadWeights; /* Weights in the standard element for Gaussian quadrature */
671
double *quadShapeFuncs; /* Shape functions evaluated at quadrature points */
672
double *quadShapeFuncDers; /* Shape function derivatives evaluated at quadrature points */
673
double x, y; /* The integration point */
674
double dxxi; /* \PartDer{x}{\xi} */
675
double dxet; /* \PartDer{x}{\eta} */
676
double dyxi; /* \PartDer{y}{\xi} */
677
double dyet; /* \PartDer{y}{\eta} */
678
double dxix; /* \PartDer{\xi}{x} */
679
double detx; /* \PartDer{\eta}{x} */
680
double dxiy; /* \PartDer{\xi}{y} */
681
double dety; /* \PartDer{\eta}{y} */
682
PetscScalar dfxi; /* \PartDer{field}{\xi} */
683
PetscScalar dfet; /* \PartDer{field}{\eta} */
684
double coords[12]; /* Coordinates of the element */
686
int i, j, k, func, p, arg;
687
#ifdef PETSC_USE_BOPT_g
693
ierr = MPI_Comm_rank(disc->comm, &rank); CHKERRQ(ierr);
694
numQuadPoints = disc->numQuadPoints;
695
quadWeights = disc->quadWeights;
696
quadShapeFuncs = disc->quadShapeFuncs;
697
quadShapeFuncDers = disc->quadShapeFuncDers;
699
/* Calculate the determinant of the inverse Jacobian of the map to the standard element */
700
if (mesh->isPeriodic == PETSC_TRUE) {
701
coords[0*2+0] = nodes[elements[elem*numCorners+0]*2+0];
702
coords[0*2+1] = nodes[elements[elem*numCorners+0]*2+1];
703
for(func = 1; func < funcs; func++) {
704
coords[func*2+0] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+func]*2+0], coords[0*2+0]);
705
coords[func*2+1] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+func]*2+1], coords[0*2+1]);
708
for(func = 0; func < funcs; func++) {
709
coords[func*2+0] = nodes[elements[elem*numCorners+func]*2+0];
710
coords[func*2+1] = nodes[elements[elem*numCorners+func]*2+1];
713
/* Check for constant jacobian here */
715
jac = PetscAbsReal((coords[2] - coords[0])*(coords[5] - coords[1]) - (coords[4] - coords[0])*(coords[3] - coords[1]));
717
PetscPrintf(PETSC_COMM_SELF, "[%d], elem: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
718
rank, elem, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
719
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
722
#ifdef PETSC_USE_BOPT_g
723
ierr = PetscOptionsHasName(PETSC_NULL, "-trace_assembly", &opt); CHKERRQ(ierr);
724
if (opt == PETSC_TRUE) {
725
PetscPrintf(PETSC_COMM_SELF, "[%d]elem: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
726
rank, elem, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
727
PetscPrintf(PETSC_COMM_SELF, " x4: %g y4: %g x5: %g y5: %g x6: %g y6: %g\n",
728
coords[6], coords[7], coords[8], coords[9], coords[10], coords[11]);
732
/* Calculate element vector entries by Gaussian quadrature */
733
for(p = 0; p < numQuadPoints; p++) {
734
/* x = \sum^{funcs}_{f=1} x_f \phi^f(p)
735
\PartDer{x}{\xi}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\xi}
736
\PartDer{x}{\eta}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\eta}
737
y = \sum^{funcs}_{f=1} y_f \phi^f(p)
738
\PartDer{y}{\xi}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\xi}
739
\PartDer{y}{\eta}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\eta}
740
u^i = \sum^{funcs}_{f=1} u^i_f \phi^f(p)
741
\PartDer{u^i}{\xi}(p) = \sum^{funcs}_{f=1} u^i_f \PartDer{\phi^f(p)}{\xi}
742
\PartDer{u^i}{\eta}(p) = \sum^{funcs}_{f=1} u^i_f \PartDer{\phi^f(p)}{\eta} */
744
dxxi = 0.0; dyxi = 0.0;
745
dxet = 0.0; dyet = 0.0;
746
for(arg = 0; arg < numArgs; arg++)
747
for(j = 0; j < comp*3; j++)
748
fieldVal[arg][j] = 0.0;
749
for(func = 0; func < funcs; func++)
751
x += coords[func*2] *quadShapeFuncs[p*funcs+func];
752
dxxi += coords[func*2] *quadShapeFuncDers[p*funcs*2+func*2];
753
dxet += coords[func*2] *quadShapeFuncDers[p*funcs*2+func*2+1];
754
y += coords[func*2+1]*quadShapeFuncs[p*funcs+func];
755
dyxi += coords[func*2+1]*quadShapeFuncDers[p*funcs*2+func*2];
756
dyet += coords[func*2+1]*quadShapeFuncDers[p*funcs*2+func*2+1];
757
for(arg = 0; arg < numArgs; arg++) {
758
for(j = 0; j < comp; j++) {
759
fieldVal[arg][j*3] += (field[arg][func*comp+j] - ALEfield[func*comp+j])*quadShapeFuncs[p*funcs+func];
760
fieldVal[arg][j*3+1] += field[arg][func*comp+j]*quadShapeFuncDers[p*funcs*2+func*2];
761
fieldVal[arg][j*3+2] += field[arg][func*comp+j]*quadShapeFuncDers[p*funcs*2+func*2+1];
765
if (mesh->isPeriodic == PETSC_TRUE) {
766
x = MeshPeriodicX(mesh, x);
767
y = MeshPeriodicY(mesh, y);
769
jac = PetscAbsReal(dxxi*dyet - dxet*dyxi);
771
PetscPrintf(PETSC_COMM_SELF, "[%d]p: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
772
rank, p, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
773
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
775
/* These are the elements of the inverse matrix */
782
/* Convert the field derivatives to old coordinates */
783
for(arg = 0; arg < numArgs; arg++) {
784
for(j = 0; j < comp; j++) {
785
dfxi = fieldVal[arg][j*3+1];
786
dfet = fieldVal[arg][j*3+2];
787
fieldVal[arg][j*3+1] = dfxi*dxix + dfet*detx;
788
fieldVal[arg][j*3+2] = dfxi*dxiy + dfet*dety;
792
ierr = (*f)(1, comp, &x, &y, PETSC_NULL, numArgs, fieldVal, funcVal, ctx); CHKERRQ(ierr);
793
#ifdef PETSC_USE_BOPT_g
794
ierr = PetscOptionsHasName(PETSC_NULL, "-trace_assembly", &opt); CHKERRQ(ierr);
795
if (opt == PETSC_TRUE) {
796
PetscPrintf(PETSC_COMM_SELF, "[%d]p: %d jac: %g", rank, p, jac);
797
for(j = 0; j < comp; j++)
798
PetscPrintf(PETSC_COMM_SELF, " func[%d]: %g", j, PetscRealPart(funcVal[j]));
799
PetscPrintf(PETSC_COMM_SELF, "\n");
803
for(i = 0, k = 0; i < funcs; i++) {
804
for(j = 0; j < comp; j++, k++) {
805
vec[k] += alpha*funcVal[j]*quadShapeFuncs[p*funcs+i]*jac*quadWeights[p];
806
#ifdef PETSC_USE_BOPT_g
807
ierr = PetscOptionsHasName(PETSC_NULL, "-trace_assembly", &opt); CHKERRQ(ierr);
808
if (opt == PETSC_TRUE) {
809
PetscPrintf(PETSC_COMM_SELF, "[%d] vec[%d]: %g\n", rank, k, PetscRealPart(vec[k]));
815
PetscLogFlops(((12 + (7*numArgs + 5)*comp)*funcs + 8 + 6*numArgs*comp) * numQuadPoints);
816
PetscFunctionReturn(0);
820
#define __FUNCT__ "Identity_Triangular_2D_Quadratic"
821
int Identity_Triangular_2D_Quadratic(Discretization disc, Discretization test, int rowSize, int colSize,
822
int globalRowStart, int globalColStart, int globalSize, double *coords,
823
PetscScalar alpha, PetscScalar *field, PetscScalar *array, void *ctx)
825
int numQuadPoints; /* Number of points used for Gaussian quadrature */
826
double *quadWeights; /* Weights in the standard element for Gaussian quadrature */
827
double *quadShapeFuncs; /* Shape functions evaluated at quadrature points */
828
double *quadTestFuncs; /* Test functions evaluated at quadrature points */
829
double *quadShapeFuncDers; /* Shape function derivatives evaluated at quadrature points */
830
double dxxi; /* \PartDer{x}{\xi} */
831
double dxet; /* \PartDer{x}{\eta} */
832
double dyxi; /* \PartDer{y}{\xi} */
833
double dyet; /* \PartDer{y}{\eta} */
834
double jac; /* |J| for map to standard element */
835
int comp; /* The number of field components */
836
int funcs; /* The number of shape functions */
840
/* Calculate element matrix entries by Gaussian quadrature */
843
numQuadPoints = disc->numQuadPoints;
844
quadWeights = disc->quadWeights;
845
quadShapeFuncs = disc->quadShapeFuncs;
846
quadTestFuncs = test->quadShapeFuncs;
847
quadShapeFuncDers = disc->quadShapeFuncDers;
848
for(p = 0; p < numQuadPoints; p++)
850
/* \PartDer{x}{\xi}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\xi}
851
\PartDer{x}{\eta}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\eta}
852
\PartDer{y}{\xi}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\xi}
853
\PartDer{y}{\eta}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\eta} */
854
dxxi = 0.0; dxet = 0.0;
855
dyxi = 0.0; dyet = 0.0;
856
for(f = 0; f < funcs; f++)
858
dxxi += coords[f*2] *quadShapeFuncDers[p*funcs*2+f*2];
859
dxet += coords[f*2] *quadShapeFuncDers[p*funcs*2+f*2+1];
860
dyxi += coords[f*2+1]*quadShapeFuncDers[p*funcs*2+f*2];
861
dyet += coords[f*2+1]*quadShapeFuncDers[p*funcs*2+f*2+1];
863
jac = PetscAbsReal(dxxi*dyet - dxet*dyxi);
864
#ifdef PETSC_USE_BOPT_g
866
PetscPrintf(PETSC_COMM_SELF, "p: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
867
p, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
868
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
872
for(i = 0; i < funcs; i++)
874
for(j = 0; j < funcs; j++)
876
for(c = 0; c < comp; c++)
878
array[(i*comp+c+globalRowStart)*globalSize + j*comp+c+globalColStart] +=
879
alpha*quadTestFuncs[p*funcs+i]*quadShapeFuncs[p*funcs+j]*jac*quadWeights[p];
880
/* PetscPrintf(PETSC_COMM_SELF, " array[%d,%d]: %g\n", i*comp+c+globalRowStart, j*comp+c+globalColStart,
881
array[(i*comp+c+globalRowStart)*globalSize + j*comp+c+globalColStart]); */
886
PetscLogFlops((8*funcs + 3 + 5*funcs*funcs*comp) * numQuadPoints);
888
PetscFunctionReturn(0);
892
#define __FUNCT__ "Laplacian_Triangular_2D_Quadratic"
893
int Laplacian_Triangular_2D_Quadratic(Discretization disc, Discretization test, int rowSize, int colSize,
894
int globalRowStart, int globalColStart, int globalSize, double *coords,
895
PetscScalar alpha, PetscScalar *field, PetscScalar *array, void *ctx)
897
int numQuadPoints; /* Number of points used for Gaussian quadrature */
898
double *quadWeights; /* Weights in the standard element for Gaussian quadrature */
899
double *quadShapeFuncDers; /* Shape function derivatives evaluated at quadrature points */
900
double *quadTestFuncDers; /* Test function derivatives evaluated at quadrature points */
901
double dxxi; /* \PartDer{x}{\xi} */
902
double dxet; /* \PartDer{x}{\eta} */
903
double dyxi; /* \PartDer{y}{\xi} */
904
double dyet; /* \PartDer{y}{\eta} */
905
double dxix; /* \PartDer{\xi}{x} */
906
double detx; /* \PartDer{\eta}{x} */
907
double dxiy; /* \PartDer{\xi}{y} */
908
double dety; /* \PartDer{\eta}{y} */
909
double dphix; /* \PartDer{\phi_i}{x} \times \PartDer{\phi_j}{x} */
910
double dphiy; /* \PartDer{\phi_i}{y} \times \PartDer{\phi_j}{y} */
911
double jac; /* |J| for map to standard element */
912
double invjac; /* |J^{-1}| for map from standard element */
913
int comp; /* The number of field components */
914
int funcs; /* The number of shape functions */
918
/* Calculate element matrix entries by Gaussian quadrature */
921
numQuadPoints = disc->numQuadPoints;
922
quadWeights = disc->quadWeights;
923
quadShapeFuncDers = disc->quadShapeFuncDers;
924
quadTestFuncDers = test->quadShapeFuncDers;
925
for(p = 0; p < numQuadPoints; p++) {
926
/* \PartDer{x}{\xi}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\xi}
927
\PartDer{x}{\eta}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\eta}
928
\PartDer{y}{\xi}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\xi}
929
\PartDer{y}{\eta}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\eta} */
930
dxxi = 0.0; dxet = 0.0;
931
dyxi = 0.0; dyet = 0.0;
932
for(f = 0; f < funcs; f++) {
933
dxxi += coords[f*2] *quadShapeFuncDers[p*funcs*2+f*2];
934
dxet += coords[f*2] *quadShapeFuncDers[p*funcs*2+f*2+1];
935
dyxi += coords[f*2+1]*quadShapeFuncDers[p*funcs*2+f*2];
936
dyet += coords[f*2+1]*quadShapeFuncDers[p*funcs*2+f*2+1];
938
jac = PetscAbsReal(dxxi*dyet - dxet*dyxi);
939
#ifdef PETSC_USE_BOPT_g
941
PetscPrintf(PETSC_COMM_SELF, "p: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
942
p, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
943
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
946
/* These are the elements of the inverse matrix */
953
for(i = 0; i < funcs; i++) {
954
for(j = 0; j < funcs; j++) {
955
dphix = (quadTestFuncDers[p*funcs*2+i*2]*dxix + quadTestFuncDers[p*funcs*2+i*2+1]*detx)*
956
(quadShapeFuncDers[p*funcs*2+j*2]*dxix + quadShapeFuncDers[p*funcs*2+j*2+1]*detx);
957
dphiy = (quadTestFuncDers[p*funcs*2+i*2]*dxiy + quadTestFuncDers[p*funcs*2+i*2+1]*dety)*
958
(quadShapeFuncDers[p*funcs*2+j*2]*dxiy + quadShapeFuncDers[p*funcs*2+j*2+1]*dety);
959
for(c = 0; c < comp; c++) {
960
array[(i*comp+c+globalRowStart)*globalSize + j*comp+c+globalColStart] +=
961
-alpha*(dphix + dphiy)*jac*quadWeights[p];
962
/* PetscPrintf(PETSC_COMM_SELF, " array[%d,%d]: %g\n", i*comp+c+globalRowStart, j*comp+c+globalColStart,
963
array[(i*comp+c+globalRowStart)*globalSize + j*comp+c+globalColStart]); */
968
PetscLogFlops((8*funcs + 8 + 19*funcs*funcs*comp) * numQuadPoints);
970
PetscFunctionReturn(0);
974
#define __FUNCT__ "Weighted_Laplacian_Triangular_2D_Quadratic"
975
int Weighted_Laplacian_Triangular_2D_Quadratic(Discretization disc, Discretization test, int rowSize, int colSize,
976
int globalRowStart, int globalColStart, int globalSize, double *coords,
977
PetscScalar alpha, PetscScalar *field, PetscScalar *array, void *ctx)
979
int numQuadPoints; /* Number of points used for Gaussian quadrature */
980
double *quadWeights; /* Weights in the standard element for Gaussian quadrature */
981
double *quadShapeFuncDers; /* Shape function derivatives evaluated at quadrature points */
982
double *quadTestFuncDers; /* Test function derivatives evaluated at quadrature points */
983
double dxxi; /* \PartDer{x}{\xi} */
984
double dxet; /* \PartDer{x}{\eta} */
985
double dyxi; /* \PartDer{y}{\xi} */
986
double dyet; /* \PartDer{y}{\eta} */
987
double dxix; /* \PartDer{\xi}{x} */
988
double detx; /* \PartDer{\eta}{x} */
989
double dxiy; /* \PartDer{\xi}{y} */
990
double dety; /* \PartDer{\eta}{y} */
991
double dphix; /* \PartDer{\phi_i}{x} \times \PartDer{\phi_j}{x} */
992
double dphiy; /* \PartDer{\phi_i}{y} \times \PartDer{\phi_j}{y} */
993
double jac; /* |J| for map to standard element */
994
double invjac; /* |J^{-1}| for map from standard element */
995
int comp; /* The number of field components */
996
int funcs; /* The number of shape functions */
999
/* Each element is weighted by its Jacobian. This is supposed to make smaller elements "stiffer". */
1001
/* Calculate element matrix entries by Gaussian quadrature */
1003
funcs = disc->funcs;
1004
numQuadPoints = disc->numQuadPoints;
1005
quadWeights = disc->quadWeights;
1006
quadShapeFuncDers = disc->quadShapeFuncDers;
1007
quadTestFuncDers = test->quadShapeFuncDers;
1008
for(p = 0; p < numQuadPoints; p++)
1010
/* \PartDer{x}{\xi}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\xi}
1011
\PartDer{x}{\eta}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\eta}
1012
\PartDer{y}{\xi}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\xi}
1013
\PartDer{y}{\eta}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\eta} */
1014
dxxi = 0.0; dxet = 0.0;
1015
dyxi = 0.0; dyet = 0.0;
1016
for(f = 0; f < funcs; f++)
1018
dxxi += coords[f*2] *quadShapeFuncDers[p*funcs*2+f*2];
1019
dxet += coords[f*2] *quadShapeFuncDers[p*funcs*2+f*2+1];
1020
dyxi += coords[f*2+1]*quadShapeFuncDers[p*funcs*2+f*2];
1021
dyet += coords[f*2+1]*quadShapeFuncDers[p*funcs*2+f*2+1];
1023
jac = PetscAbsReal(dxxi*dyet - dxet*dyxi);
1024
#ifdef PETSC_USE_BOPT_g
1025
if (jac < 1.0e-14) {
1026
PetscPrintf(PETSC_COMM_SELF, "p: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
1027
p, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
1028
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
1031
/* These are the elements of the inverse matrix */
1034
dxiy = -dxet*invjac;
1035
detx = -dyxi*invjac;
1038
for(i = 0; i < funcs; i++)
1040
for(j = 0; j < funcs; j++)
1042
dphix = (quadTestFuncDers[p*funcs*2+i*2]*dxix + quadTestFuncDers[p*funcs*2+i*2+1]*detx)*
1043
(quadShapeFuncDers[p*funcs*2+j*2]*dxix + quadShapeFuncDers[p*funcs*2+j*2+1]*detx);
1044
dphiy = (quadTestFuncDers[p*funcs*2+i*2]*dxiy + quadTestFuncDers[p*funcs*2+i*2+1]*dety)*
1045
(quadShapeFuncDers[p*funcs*2+j*2]*dxiy + quadShapeFuncDers[p*funcs*2+j*2+1]*dety);
1046
for(c = 0; c < comp; c++)
1048
array[(i*comp+c+globalRowStart)*globalSize + j*comp+c+globalColStart] +=
1049
-alpha*(dphix + dphiy)*quadWeights[p];
1054
PetscLogFlops((8*funcs + 8 + 18*funcs*funcs*comp) * numQuadPoints);
1056
PetscFunctionReturn(0);
1060
#define __FUNCT__ "Divergence_Triangular_2D_Quadratic"
1061
int Divergence_Triangular_2D_Quadratic(Discretization disc, Discretization test, int rowSize, int colSize,
1062
int globalRowStart, int globalColStart, int globalSize, double *coords,
1063
PetscScalar alpha, PetscScalar *field, PetscScalar *array, void *ctx)
1065
/* We are using the convention that
1067
\nabla \matrix{v_1 \cr v_2 \cr \vdots \cr v_n} =
1068
\matrix{v^{(1)}_1 \cr \vdots \cr v^{(d)}_1 \cr v^{(1)}_2 \cr \vdots \cr v^{(d)}_n}
1072
\nabla \cdot \matrix{v^{(1)}_1 \cr \vdots \cr v^{(d)}_1 \cr v^{(1)}_2 \cr \vdots \cr v^{(d)}_n} =
1073
\matrix{v_1 \cr v_2 \cr \vdots \cr v_n}
1075
where $d$ is the number of space dimensions. This agrees with the convention which allows
1076
$\Delta \matrix{u_1 \cr u_2} = 0$ to denote a set of scalar equations This also requires that
1077
the dimension of a vector must be divisible by the space dimension in order to be acted upon by
1078
the divergence operator */
1079
int numQuadPoints; /* Number of points used for Gaussian quadrature */
1080
double *quadWeights; /* Weights in the standard element for Gaussian quadrature */
1081
double *quadTestFuncs; /* Test functions evaluated at quadrature points */
1082
double *quadShapeFuncDers; /* Shape function derivatives evaluated at quadrature points */
1083
double dxxi; /* \PartDer{x}{\xi} */
1084
double dxet; /* \PartDer{x}{\eta} */
1085
double dyxi; /* \PartDer{y}{\xi} */
1086
double dyet; /* \PartDer{y}{\eta} */
1087
double dxix; /* \PartDer{\xi}{x} */
1088
double detx; /* \PartDer{\eta}{x} */
1089
double dxiy; /* \PartDer{\xi}{y} */
1090
double dety; /* \PartDer{\eta}{y} */
1091
double dphix; /* \PartDer{\phi_i}{x} \times \PartDer{\phi_j}{x} */
1092
double dphiy; /* \PartDer{\phi_i}{y} \times \PartDer{\phi_j}{y} */
1093
double jac; /* |J| for map to standard element */
1094
double invjac; /* |J^{-1}| for map from standard element */
1095
int comp; /* The number of field components */
1096
int tcomp; /* The number of field components for the test field */
1097
int funcs; /* The number of shape functions */
1098
int tfuncs; /* The number of test functions */
1099
int i, j, c, tc, f, p;
1102
/* Calculate element matrix entries by Gaussian quadrature */
1105
funcs = disc->funcs;
1106
tfuncs = test->funcs;
1107
numQuadPoints = disc->numQuadPoints;
1108
quadWeights = disc->quadWeights;
1109
quadTestFuncs = test->quadShapeFuncs;
1110
quadShapeFuncDers = disc->quadShapeFuncDers;
1111
for(p = 0; p < numQuadPoints; p++) {
1112
/* \PartDer{x}{\xi}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\xi}
1113
\PartDer{x}{\eta}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\eta}
1114
\PartDer{y}{\xi}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\xi}
1115
\PartDer{y}{\eta}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\eta} */
1116
dxxi = 0.0; dxet = 0.0;
1117
dyxi = 0.0; dyet = 0.0;
1118
for(f = 0; f < funcs; f++)
1120
dxxi += coords[f*2] *quadShapeFuncDers[p*funcs*2+f*2];
1121
dxet += coords[f*2] *quadShapeFuncDers[p*funcs*2+f*2+1];
1122
dyxi += coords[f*2+1]*quadShapeFuncDers[p*funcs*2+f*2];
1123
dyet += coords[f*2+1]*quadShapeFuncDers[p*funcs*2+f*2+1];
1125
jac = PetscAbsReal(dxxi*dyet - dxet*dyxi);
1126
#ifdef PETSC_USE_BOPT_g
1127
if (jac < 1.0e-14) {
1128
PetscPrintf(PETSC_COMM_SELF, "p: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
1129
p, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
1130
SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
1133
/* PetscPrintf(PETSC_COMM_SELF, "p: %d x1: %g y1: %g x2: %g y2: %g x3: %g y3: %g\n",
1134
p, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]); */
1135
/* These are the elements of the inverse matrix */
1138
dxiy = -dxet*invjac;
1139
detx = -dyxi*invjac;
1142
/* The rows are test functions */
1143
for(i = 0; i < tfuncs; i++)
1145
for(tc = 0; tc < tcomp; tc++)
1147
/* The columns are shape functions */
1148
for(j = 0; j < funcs; j++)
1150
dphix = quadShapeFuncDers[p*funcs*2+j*2]*dxix + quadShapeFuncDers[p*funcs*2+j*2+1]*detx;
1151
dphiy = quadShapeFuncDers[p*funcs*2+j*2]*dxiy + quadShapeFuncDers[p*funcs*2+j*2+1]*dety;
1152
/* We divide by the number of space dimensions */
1153
for(c = 0; c < comp/2; c++)
1155
array[(i*tcomp+tc+globalRowStart)*globalSize + j*comp+c*2+globalColStart] +=
1156
alpha*dphix*quadTestFuncs[p*tfuncs+i]*jac*quadWeights[p];
1157
/* PetscPrintf(PETSC_COMM_SELF, " array[%d,%d]: %g\n", i*tcomp+tc+globalRowStart, j*comp+c*2+globalColStart,
1158
array[(i*tcomp+tc+globalRowStart)*globalSize + j*comp+c*2+globalColStart]); */
1159
array[(i*tcomp+tc+globalRowStart)*globalSize + j*comp+c*2+1+globalColStart] +=
1160
alpha*dphiy*quadTestFuncs[p*tfuncs+i]*jac*quadWeights[p];
1161
/* PetscPrintf(PETSC_COMM_SELF, " array[%d,%d]: %g\n", i*tcomp+tc+globalRowStart, j*comp+c*2+1+globalColStart,
1162
array[(i*tcomp+tc+globalRowStart)*globalSize + j*comp+c*2+1+globalColStart]); */
1168
PetscLogFlops((8*funcs + 8 + 8*tfuncs*tcomp*funcs*comp) * numQuadPoints);
1170
PetscFunctionReturn(0);
1174
#define __FUNCT__ "DiscInterpolateField_Triangular_2D_Quadratic"
1175
int DiscInterpolateField_Triangular_2D_Quadratic(Discretization disc, Mesh oldMesh, int elem, double x, double y, double z,
1176
PetscScalar *oldFieldVal, PetscScalar *newFieldVal, InterpolationType type)
1178
Mesh_Triangular *tri = (Mesh_Triangular *) oldMesh->data;
1179
int numCorners = oldMesh->numCorners;
1180
int *elements = tri->faces;
1181
int *neighbors = tri->neighbors;
1182
double *nodes = tri->nodes;
1183
double coords[24]; /* Coordinates of our "big element" */
1184
double xi, eta; /* Canonical coordinates of the point */
1185
double x21, x31, y21, y31, jac, invjac, dx, dy, dxix, dxiy, detx, dety, xiOld, etaOld;
1186
int comp = disc->comp;
1187
int neighbor, corner, nelem, node, c;
1191
/* No scheme in place for boundary elements */
1192
for(neighbor = 0; neighbor < 3; neighbor++)
1193
if (neighbors[elem*3+neighbor] < 0)
1194
type = INTERPOLATION_LOCAL;
1198
case INTERPOLATION_LOCAL:
1199
if (oldMesh->isPeriodic == PETSC_TRUE) {
1200
coords[0*2+0] = MeshPeriodicRelativeX(oldMesh, nodes[elements[elem*numCorners+0]*2+0], x);
1201
coords[0*2+1] = MeshPeriodicRelativeY(oldMesh, nodes[elements[elem*numCorners+0]*2+1], y);
1202
coords[1*2+0] = MeshPeriodicRelativeX(oldMesh, nodes[elements[elem*numCorners+1]*2+0], x);
1203
coords[1*2+1] = MeshPeriodicRelativeY(oldMesh, nodes[elements[elem*numCorners+1]*2+1], y);
1204
coords[2*2+0] = MeshPeriodicRelativeX(oldMesh, nodes[elements[elem*numCorners+2]*2+0], x);
1205
coords[2*2+1] = MeshPeriodicRelativeY(oldMesh, nodes[elements[elem*numCorners+2]*2+1], y);
1206
coords[3*2+0] = MeshPeriodicRelativeX(oldMesh, nodes[elements[elem*numCorners+3]*2+0], x);
1207
coords[3*2+1] = MeshPeriodicRelativeY(oldMesh, nodes[elements[elem*numCorners+3]*2+1], y);
1208
coords[4*2+0] = MeshPeriodicRelativeX(oldMesh, nodes[elements[elem*numCorners+4]*2+0], x);
1209
coords[4*2+1] = MeshPeriodicRelativeY(oldMesh, nodes[elements[elem*numCorners+4]*2+1], y);
1210
coords[5*2+0] = MeshPeriodicRelativeX(oldMesh, nodes[elements[elem*numCorners+5]*2+0], x);
1211
coords[5*2+1] = MeshPeriodicRelativeY(oldMesh, nodes[elements[elem*numCorners+5]*2+1], y);
1213
coords[0*2+0] = nodes[elements[elem*numCorners+0]*2+0];
1214
coords[0*2+1] = nodes[elements[elem*numCorners+0]*2+1];
1215
coords[1*2+0] = nodes[elements[elem*numCorners+1]*2+0];
1216
coords[1*2+1] = nodes[elements[elem*numCorners+1]*2+1];
1217
coords[2*2+0] = nodes[elements[elem*numCorners+2]*2+0];
1218
coords[2*2+1] = nodes[elements[elem*numCorners+2]*2+1];
1219
coords[3*2+0] = nodes[elements[elem*numCorners+3]*2+0];
1220
coords[3*2+1] = nodes[elements[elem*numCorners+3]*2+1];
1221
coords[4*2+0] = nodes[elements[elem*numCorners+4]*2+0];
1222
coords[4*2+1] = nodes[elements[elem*numCorners+4]*2+1];
1223
coords[5*2+0] = nodes[elements[elem*numCorners+5]*2+0];
1224
coords[5*2+1] = nodes[elements[elem*numCorners+5]*2+1];
1226
/* Get the (\xi,\eta) coordinates of the point */
1227
ierr = DiscTransformCoords_Triangular_2D_Quadratic(x, y, coords, &xi, &eta); CHKERRQ(ierr);
1228
if ((xi < -1.0e-02) || (eta < -1.0e-02) || (xi > 1.01) || (eta > 1.01)) {
1231
/* Use linear approximation */
1232
x21 = coords[1*2+0] - coords[0*2+0];
1233
x31 = coords[2*2+0] - coords[0*2+0];
1234
y21 = coords[1*2+1] - coords[0*2+1];
1235
y31 = coords[2*2+1] - coords[0*2+1];
1236
dx = x - coords[0*2+0];
1237
dy = y - coords[0*2+1];
1238
jac = PetscAbsReal(x21*y31 - x31*y21);
1245
xi = dxix*dx + dxiy*dy;
1246
eta = detx*dx + dety*dy;
1247
PetscPrintf(PETSC_COMM_SELF, "elem: %d x: %g y: %g xiOld: %g etaOld: %g xi: %g eta: %g\n", elem, x, y, xiOld, etaOld, xi, eta);
1249
for(c = 0; c < comp; c++) {
1250
newFieldVal[c] = oldFieldVal[0*comp+c]*(1.0 - xi - eta)*(1.0 - 2.0*xi - 2.0*eta) +
1251
oldFieldVal[1*comp+c]*xi *(2.0*xi - 1.0) +
1252
oldFieldVal[2*comp+c]*eta*(2.0*eta - 1.0) +
1253
oldFieldVal[3*comp+c]*4.0*xi*eta +
1254
oldFieldVal[4*comp+c]*4.0*eta*(1.0 - xi - eta) +
1255
oldFieldVal[5*comp+c]*4.0*xi *(1.0 - xi - eta);
1257
PetscLogFlops(34*comp);
1259
case INTERPOLATION_HALO:
1260
/* Here is our "big element" where numbers in parantheses represent
1261
the numbering on the old little element:
1282
We search for the neighbor node by looking for the vertex not a member of the original element.
1284
for(neighbor = 0; neighbor < 3; neighbor++)
1286
nelem = neighbors[elem*3+neighbor];
1287
for(corner = 0; corner < 3; corner++)
1289
node = elements[nelem*numCorners+corner];
1290
if ((node != elements[elem*numCorners+((neighbor+1)%3)]) && (node != elements[elem*numCorners+((neighbor+2)%3)]))
1292
/* The neighboring elements give the vertices */
1293
coords[neighbor*2] = nodes[node*2];
1294
coords[neighbor*2+1] = nodes[node*2+1];
1299
/* Element vertices form midnodes */
1300
coords[3*2] = nodes[elements[elem*numCorners]*2];
1301
coords[3*2+1] = nodes[elements[elem*numCorners]*2+1];
1302
coords[4*2] = nodes[elements[elem*numCorners+1]*2];
1303
coords[4*2+1] = nodes[elements[elem*numCorners+1]*2+1];
1304
coords[5*2] = nodes[elements[elem*numCorners+2]*2];
1305
coords[5*2+1] = nodes[elements[elem*numCorners+2]*2+1];
1306
/* Treat 4 triangles as one big element with quadratic shape functions */
1307
SETERRQ(PETSC_ERR_SUP, "Unsupported interpolation type");
1309
SETERRQ(PETSC_ERR_ARG_WRONG, "Unknown interpolation type");
1312
PetscFunctionReturn(0);
1316
#define __FUNCT__ "DiscInterpolateElementVec_Triangular_2D_Quadratic"
1317
int DiscInterpolateElementVec_Triangular_2D_Quadratic(Discretization disc, ElementVec vec, Discretization newDisc, ElementVec newVec)
1319
int comp = disc->comp;
1320
int size = disc->size;
1321
PetscScalar *array, *newArray;
1322
PetscTruth islin, isquad;
1327
ierr = ElementVecGetArray(vec, &array); CHKERRQ(ierr);
1328
ierr = ElementVecGetArray(newVec, &newArray); CHKERRQ(ierr);
1329
ierr = PetscTypeCompare((PetscObject) newDisc, DISCRETIZATION_TRIANGULAR_2D_LINEAR, &islin); CHKERRQ(ierr);
1330
ierr = PetscTypeCompare((PetscObject) newDisc, DISCRETIZATION_TRIANGULAR_2D_QUADRATIC, &isquad); CHKERRQ(ierr);
1331
if (isquad == PETSC_TRUE) {
1332
ierr = PetscMemcpy(newArray, array, size * sizeof(PetscScalar)); CHKERRQ(ierr);
1333
} else if (islin == PETSC_TRUE) {
1334
for(f = 0; f < newDisc->funcs; f++) {
1335
for(c = 0; c < comp; c++) {
1336
newArray[f*comp+c] = array[f*comp+c];
1340
SETERRQ(PETSC_ERR_SUP, "Discretization not supported");
1342
ierr = ElementVecRestoreArray(vec, &array); CHKERRQ(ierr);
1343
ierr = ElementVecRestoreArray(newVec, &newArray); CHKERRQ(ierr);
1344
PetscFunctionReturn(0);
1348
#define __FUNCT__ "DiscSetupQuadrature_Triangular_2D_Quadratic"
1350
DiscSetupQuadrature_Triangular_2D_Quadratic - Setup Gaussian quadrature with a 7 point integration rule
1353
. disc - The Discretization
1355
int DiscSetupQuadrature_Triangular_2D_Quadratic(Discretization disc) {
1356
int dim = disc->dim;
1357
int funcs = disc->funcs;
1363
disc->numQuadPoints = 7;
1364
ierr = PetscMalloc(disc->numQuadPoints*dim * sizeof(double), &disc->quadPoints); CHKERRQ(ierr);
1365
ierr = PetscMalloc(disc->numQuadPoints * sizeof(double), &disc->quadWeights); CHKERRQ(ierr);
1366
ierr = PetscMalloc(disc->numQuadPoints*funcs * sizeof(double), &disc->quadShapeFuncs); CHKERRQ(ierr);
1367
ierr = PetscMalloc(disc->numQuadPoints*funcs*dim * sizeof(double), &disc->quadShapeFuncDers); CHKERRQ(ierr);
1368
PetscLogObjectMemory(disc, (disc->numQuadPoints*(funcs*(dim+1) + dim+1)) * sizeof(double));
1369
disc->quadPoints[0] = 1.0/3.0;
1370
disc->quadPoints[1] = disc->quadPoints[0];
1371
disc->quadWeights[0] = 0.11250000000000;
1372
disc->quadPoints[2] = 0.797426985353087;
1373
disc->quadPoints[3] = 0.101286507323456;
1374
disc->quadWeights[1] = 0.0629695902724135;
1375
disc->quadPoints[4] = disc->quadPoints[3];
1376
disc->quadPoints[5] = disc->quadPoints[2];
1377
disc->quadWeights[2] = disc->quadWeights[1];
1378
disc->quadPoints[6] = disc->quadPoints[4];
1379
disc->quadPoints[7] = disc->quadPoints[3];
1380
disc->quadWeights[3] = disc->quadWeights[1];
1381
disc->quadPoints[8] = 0.470142064105115;
1382
disc->quadPoints[9] = 0.059715871789770;
1383
disc->quadWeights[4] = 0.066197076394253;
1384
disc->quadPoints[10] = disc->quadPoints[8];
1385
disc->quadPoints[11] = disc->quadPoints[8];
1386
disc->quadWeights[5] = disc->quadWeights[4];
1387
disc->quadPoints[12] = disc->quadPoints[9];
1388
disc->quadPoints[13] = disc->quadPoints[8];
1389
disc->quadWeights[6] = disc->quadWeights[4];
1390
for(p = 0; p < disc->numQuadPoints; p++) {
1391
xi = disc->quadPoints[p*2];
1392
eta = disc->quadPoints[p*2+1];
1393
/* \phi^0: 1 - 3 (\xi + \eta) + 2 (\xi + \eta)^2 */
1394
disc->quadShapeFuncs[p*funcs] = 1.0 - 3.0*(xi + eta) + 2.0*(xi + eta)*(xi + eta);
1395
disc->quadShapeFuncDers[p*funcs*2+0*2] = -3.0 + 4.0*(xi + eta);
1396
disc->quadShapeFuncDers[p*funcs*2+0*2+1] = -3.0 + 4.0*(xi + eta);
1397
/* \phi^1: \xi (2\xi - 1) */
1398
disc->quadShapeFuncs[p*funcs+1] = xi*(2.0*xi - 1.0);
1399
disc->quadShapeFuncDers[p*funcs*2+1*2] = 4.0*xi - 1.0;
1400
disc->quadShapeFuncDers[p*funcs*2+1*2+1] = 0.0;
1401
/* \phi^2: \eta (2\eta - 1) */
1402
disc->quadShapeFuncs[p*funcs+2] = eta*(2.0*eta - 1.0);
1403
disc->quadShapeFuncDers[p*funcs*2+2*2] = 0.0;
1404
disc->quadShapeFuncDers[p*funcs*2+2*2+1] = 4.0*eta - 1.0;
1405
/* \phi^3: 4 \xi \eta */
1406
disc->quadShapeFuncs[p*funcs+3] = 4.0*xi*eta;
1407
disc->quadShapeFuncDers[p*funcs*2+3*2] = 4.0*eta;
1408
disc->quadShapeFuncDers[p*funcs*2+3*2+1] = 4.0*xi;
1409
/* \phi^4: 4 \eta (1 - \xi - \eta) */
1410
disc->quadShapeFuncs[p*funcs+4] = 4.0*eta*(1.0 - xi - eta);
1411
disc->quadShapeFuncDers[p*funcs*2+4*2] = -4.0*eta;
1412
disc->quadShapeFuncDers[p*funcs*2+4*2+1] = -8.0*eta + 4.0*(1.0 - xi);
1413
/* \phi^5: 4 \xi (1 - \xi - \eta) */
1414
disc->quadShapeFuncs[p*funcs+5] = 4.0*xi*(1.0 - xi - eta);
1415
disc->quadShapeFuncDers[p*funcs*2+5*2] = -8.0*xi + 4.0*(1.0 - eta);
1416
disc->quadShapeFuncDers[p*funcs*2+5*2+1] = -4.0*xi;
1418
PetscFunctionReturn(0);
1422
#define __FUNCT__ "DiscSetupOperators_Triangular_2D_Quadratic"
1424
DiscSetupOperators_Triangular_2D_Quadratic - Setup the default operators
1427
. disc - The Discretization
1429
int DiscSetupOperators_Triangular_2D_Quadratic(Discretization disc) {
1434
/* The Identity operator I -- the matrix is symmetric */
1435
ierr = DiscretizationRegisterOperator(disc, Identity_Triangular_2D_Quadratic, &newOp); CHKERRQ(ierr);
1436
if (newOp != IDENTITY) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE, "Default operator %d not setup correctly", IDENTITY);
1437
/* The Laplacian operator \Delta -- the matrix is symmetric */
1438
ierr = DiscretizationRegisterOperator(disc, Laplacian_Triangular_2D_Quadratic, &newOp); CHKERRQ(ierr);
1439
if (newOp != LAPLACIAN) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE, "Default operator %d not setup correctly", LAPLACIAN);
1440
/* The Gradient operator \nabla -- the matrix is rectangular */
1441
ierr = DiscretizationRegisterOperator(disc, PETSC_NULL, &newOp); CHKERRQ(ierr);
1442
if (newOp != GRADIENT) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE, "Default operator %d not setup correctly", GRADIENT);
1443
/* The Divergence operator \nabla\cdot -- the matrix is rectangular */
1444
ierr = DiscretizationRegisterOperator(disc, Divergence_Triangular_2D_Quadratic, &newOp); CHKERRQ(ierr);
1445
if (newOp != DIVERGENCE) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE, "Default operator %d not setup correctly", DIVERGENCE);
1446
/* The weighted Laplacian operator -- the matrix is symmetric */
1447
ierr = DiscretizationRegisterOperator(disc, Weighted_Laplacian_Triangular_2D_Quadratic, &newOp); CHKERRQ(ierr);
1448
if (newOp != WEIGHTED_LAP) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE, "Default operator %d not setup correctly", WEIGHTED_LAP);
1449
PetscFunctionReturn(0);
1452
static struct _DiscretizationOps DOps = {PETSC_NULL/* DiscretizationSetup */,
1453
DiscSetupOperators_Triangular_2D_Quadratic,
1454
PETSC_NULL/* DiscretizationSetFromOptions */,
1455
DiscView_Triangular_2D_Quadratic,
1456
DiscDestroy_Triangular_2D_Quadratic,
1457
DiscEvaluateFunctionGalerkin_Triangular_2D_Quadratic,
1458
DiscEvaluateOperatorGalerkin_Triangular_2D_Quadratic,
1459
DiscEvaluateALEOperatorGalerkin_Triangular_2D_Quadratic,
1460
DiscEvaluateNonlinearOperatorGalerkin_Triangular_2D_Quadratic,
1461
DiscEvaluateNonlinearALEOperatorGalerkin_Triangular_2D_Quadratic,
1462
DiscInterpolateField_Triangular_2D_Quadratic,
1463
DiscInterpolateElementVec_Triangular_2D_Quadratic};
1467
#define __FUNCT__ "DiscCreate_Triangular_2D_Quadratic"
1468
int DiscCreate_Triangular_2D_Quadratic(Discretization disc) {
1473
if (disc->comp <= 0) {
1474
SETERRQ(PETSC_ERR_ARG_WRONG, "Discretization must have at least 1 component. Call DiscretizationSetNumComponents() to set this.");
1476
ierr = PetscMemcpy(disc->ops, &DOps, sizeof(struct _DiscretizationOps)); CHKERRQ(ierr);
1479
disc->size = disc->funcs*disc->comp;
1481
ierr = DiscretizationSetupDefaultOperators(disc); CHKERRQ(ierr);
1482
ierr = DiscSetupQuadrature_Triangular_2D_Quadratic(disc); CHKERRQ(ierr);
1484
ierr = DiscretizationCreate(disc->comm, &disc->bdDisc); CHKERRQ(ierr);
1485
ierr = DiscretizationSetNumComponents(disc->bdDisc, disc->comp); CHKERRQ(ierr);
1486
ierr = DiscretizationSetType(disc->bdDisc, BD_DISCRETIZATION_TRIANGULAR_2D_QUADRATIC); CHKERRQ(ierr);
1489
ierr = PetscMalloc(disc->comp * sizeof(PetscScalar), &disc->funcVal); CHKERRQ(ierr);
1490
ierr = PetscMalloc(2 * sizeof(PetscScalar *), &disc->fieldVal); CHKERRQ(ierr);
1491
for(arg = 0; arg < 2; arg++) {
1492
ierr = PetscMalloc(disc->comp*(disc->dim+1) * sizeof(PetscScalar), &disc->fieldVal[arg]); CHKERRQ(ierr);
1494
PetscFunctionReturn(0);
added
removed
src/grid/discretization/impls/triangular/2d/quadratic/quadratic.c
