 v0.14.0
Searching...
No Matches
l2.c
Go to the documentation of this file.
1/** \file l2.c
2
3 Based on Hierarchic Finite Element Bases on Unstructured Tetrahedral
4 Meshes, by Mark Ainsworth and Joe Coyle
5 Shape functions for MBTRI and H1 approximation
6
7*/
8
9#include <petscsys.h>
10#include <cblas.h>
11
12#include <definitions.h>
13#include <fem_tools.h>
14#include <base_functions.h>
15#include <h1_hdiv_hcurl_l2.h>
16
17static PetscErrorCode ierr;
18
20 int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM,
21 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
22 double *L, double *diffL,
23 const int dim)) {
25
26 int P = NBFACETRI_L2(p);
27 if (P == 0)
29 double diff_ksiL01, diff_ksiL20;
30 int dd = 0;
31 for (; dd < 2; dd++) {
32 diff_ksiL01[dd] = (diffN[1 * 2 + dd] - diffN[0 * 2 + dd]);
33 diff_ksiL20[dd] = (diffN[2 * 2 + dd] - diffN[0 * 2 + dd]);
34 }
35 int ii = 0;
36 for (; ii != GDIM; ++ii) {
37 int node_shift = ii * 3;
38 double ksiL01 = N[node_shift + 1] - N[node_shift + 0];
39 double ksiL20 = N[node_shift + 2] - N[node_shift + 0];
40 double L01[p + 1], L20[p + 1];
41 double diffL01[2 * (p + 1)], diffL20[2 * (p + 1)];
42 ierr = base_polynomials(p, ksiL01, diff_ksiL01, L01, diffL01, 2);
43 CHKERRQ(ierr);
44 ierr = base_polynomials(p, ksiL20, diff_ksiL20, L20, diffL20, 2);
45 CHKERRQ(ierr);
46 int shift = ii * P;
47 int jj = 0;
48 int oo = 0;
49 for (; oo <= p; oo++) {
50 int pp0 = 0;
51 for (; pp0 <= oo; pp0++) {
52 int pp1 = oo - pp0;
53 if (pp1 >= 0) {
54 if (L2N != NULL) {
55 L2N[shift + jj] = L01[pp0] * L20[pp1];
56 }
57 if (diff_L2N != NULL) {
58 int dd = 0;
59 for (; dd < 2; dd++) {
60 diff_L2N[2 * shift + 2 * jj + dd] =
61 diffL01[dd * (p + 1) + pp0] * L20[pp1] +
62 L01[pp0] * diffL20[dd * (p + 1) + pp1];
63 }
64 }
65 jj++;
66 }
67 }
68 }
69 if (jj != P)
70 SETERRQ1(PETSC_COMM_SELF, 1, "wrong order %d", jj);
71 }
73}
75 int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM,
76 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
77 double *L, double *diffL,
78 const int dim)) {
80
81 int P = NBVOLUMETET_L2(p);
82 if (P == 0)
84 double diff_ksiL0, diff_ksiL1, diff_ksiL2;
85 int dd = 0;
86 if (diffN != NULL) {
87 for (; dd < 3; dd++) {
88 diff_ksiL0[dd] = (diffN[1 * 3 + dd] - diffN[0 * 3 + dd]);
89 diff_ksiL1[dd] = (diffN[2 * 3 + dd] - diffN[0 * 3 + dd]);
90 diff_ksiL2[dd] = (diffN[3 * 3 + dd] - diffN[0 * 3 + dd]);
91 }
92 }
93 int ii = 0;
94 for (; ii != GDIM; ++ii) {
95 int node_shift = ii * 4;
96 double ksiL0 = N[node_shift + 1] - N[node_shift + 0];
97 double ksiL1 = N[node_shift + 2] - N[node_shift + 0];
98 double ksiL2 = N[node_shift + 3] - N[node_shift + 0];
99 double L0[p + 1], L1[p + 1], L2[p + 1];
100 double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
101 if (diffN != NULL) {
102 ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
103 CHKERRQ(ierr);
104 ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
105 CHKERRQ(ierr);
106 ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
107 CHKERRQ(ierr);
108 } else {
109 ierr = base_polynomials(p, ksiL0, NULL, L0, NULL, 3);
110 CHKERRQ(ierr);
111 ierr = base_polynomials(p, ksiL1, NULL, L1, NULL, 3);
112 CHKERRQ(ierr);
113 ierr = base_polynomials(p, ksiL2, NULL, L2, NULL, 3);
114 CHKERRQ(ierr);
115 }
116 int shift = ii * P;
117 int jj = 0;
118 int oo = 0;
119 for (; oo <= p; oo++) {
120 int pp0 = 0;
121 for (; pp0 <= oo; pp0++) {
122 int pp1 = 0;
123 for (; (pp0 + pp1) <= oo; pp1++) {
124 int pp2 = oo - pp0 - pp1;
125 if (pp2 >= 0) {
126 if (L2N != NULL) {
127 L2N[shift + jj] = L0[pp0] * L1[pp1] * L2[pp2];
128 }
129 if (diff_L2N != NULL) {
130 int dd = 0;
131 for (; dd < 3; dd++) {
132 diff_L2N[3 * shift + 3 * jj + dd] =
133 diffL0[dd * (p + 1) + pp0] * L1[pp1] * L2[pp2] +
134 L0[pp0] * diffL1[dd * (p + 1) + pp1] * L2[pp2] +
135 L0[pp0] * L1[pp1] * diffL2[dd * (p + 1) + pp2];
136 }
137 }
138 jj++;
139 }
140 }
141 }
142 }
143 if (jj != P)
144 SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P);
145 }
147}
148;
static Index< 'L', 3 > L
static Index< 'p', 3 > p
useful compiler directives and definitions
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const int dim
Loose implementation of some useful functions.
Functions to approximate hierarchical spaces.
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
static PetscErrorCode ierr
Definition: l2.c:17
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on tetrahedron for L2 space.
Definition: l2.c:74
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTRI(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on triangle for L2 space.
Definition: l2.c:19
const int N
Definition: speed_test.cpp:3