v0.9.0
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 /**
10  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
13  * License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
17  */
18 
19 #include <petscsys.h>
20 #include <cblas.h>
21 
22 #include <definitions.h>
23 #include <fem_tools.h>
24 #include <base_functions.h>
25 #include <h1_hdiv_hcurl_l2.h>
26 
27 static PetscErrorCode ierr;
28 
30  int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM,
31  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
32  double *L, double *diffL,
33  const int dim)) {
35 
36  int P = NBFACETRI_L2(p);
37  if (P == 0)
39  double diff_ksiL01[2], diff_ksiL20[2];
40  int dd = 0;
41  for (; dd < 2; dd++) {
42  diff_ksiL01[dd] = (diffN[1 * 2 + dd] - diffN[0 * 2 + dd]);
43  diff_ksiL20[dd] = (diffN[2 * 2 + dd] - diffN[0 * 2 + dd]);
44  }
45  int ii = 0;
46  for (; ii != GDIM; ++ii) {
47  int node_shift = ii * 3;
48  double ksiL01 = N[node_shift + 1] - N[node_shift + 0];
49  double ksiL20 = N[node_shift + 2] - N[node_shift + 0];
50  double L01[p + 1], L20[p + 1];
51  double diffL01[2 * (p + 1)], diffL20[2 * (p + 1)];
52  ierr = base_polynomials(p, ksiL01, diff_ksiL01, L01, diffL01, 2);
53  CHKERRQ(ierr);
54  ierr = base_polynomials(p, ksiL20, diff_ksiL20, L20, diffL20, 2);
55  CHKERRQ(ierr);
56  int shift = ii * P;
57  int jj = 0;
58  int oo = 0;
59  for (; oo <= p; oo++) {
60  int pp0 = 0;
61  for (; pp0 <= oo; pp0++) {
62  int pp1 = oo - pp0;
63  if (pp1 >= 0) {
64  if (L2N != NULL) {
65  L2N[shift + jj] = L01[pp0] * L20[pp1];
66  }
67  if (diff_L2N != NULL) {
68  int dd = 0;
69  for (; dd < 2; dd++) {
70  diff_L2N[2 * shift + 2 * jj + dd] =
71  diffL01[dd * (p + 1) + pp0] * L20[pp1] +
72  L01[pp0] * diffL20[dd * (p + 1) + pp1];
73  }
74  }
75  jj++;
76  }
77  }
78  }
79  if (jj != P)
80  SETERRQ1(PETSC_COMM_SELF, 1, "wrong order %d", jj);
81  }
83 }
85  int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM,
86  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
87  double *L, double *diffL,
88  const int dim)) {
90 
91  int P = NBVOLUMETET_L2(p);
92  if (P == 0)
94  double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
95  int dd = 0;
96  if (diffN != NULL) {
97  for (; dd < 3; dd++) {
98  diff_ksiL0[dd] = (diffN[1 * 3 + dd] - diffN[0 * 3 + dd]);
99  diff_ksiL1[dd] = (diffN[2 * 3 + dd] - diffN[0 * 3 + dd]);
100  diff_ksiL2[dd] = (diffN[3 * 3 + dd] - diffN[0 * 3 + dd]);
101  }
102  }
103  int ii = 0;
104  for (; ii != GDIM; ++ii) {
105  int node_shift = ii * 4;
106  double ksiL0 = N[node_shift + 1] - N[node_shift + 0];
107  double ksiL1 = N[node_shift + 2] - N[node_shift + 0];
108  double ksiL2 = N[node_shift + 3] - N[node_shift + 0];
109  double L0[p + 1], L1[p + 1], L2[p + 1];
110  double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
111  if (diffN != NULL) {
112  ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
113  CHKERRQ(ierr);
114  ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
115  CHKERRQ(ierr);
116  ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
117  CHKERRQ(ierr);
118  } else {
119  ierr = base_polynomials(p, ksiL0, NULL, L0, NULL, 3);
120  CHKERRQ(ierr);
121  ierr = base_polynomials(p, ksiL1, NULL, L1, NULL, 3);
122  CHKERRQ(ierr);
123  ierr = base_polynomials(p, ksiL2, NULL, L2, NULL, 3);
124  CHKERRQ(ierr);
125  }
126  int shift = ii * P;
127  int jj = 0;
128  int oo = 0;
129  for (; oo <= p; oo++) {
130  int pp0 = 0;
131  for (; pp0 <= oo; pp0++) {
132  int pp1 = 0;
133  for (; (pp0 + pp1) <= oo; pp1++) {
134  int pp2 = oo - pp0 - pp1;
135  if (pp2 >= 0) {
136  if (L2N != NULL) {
137  L2N[shift + jj] = L0[pp0] * L1[pp1] * L2[pp2];
138  }
139  if (diff_L2N != NULL) {
140  int dd = 0;
141  for (; dd < 3; dd++) {
142  diff_L2N[3 * shift + 3 * jj + dd] =
143  diffL0[dd * (p + 1) + pp0] * L1[pp1] * L2[pp2] +
144  L0[pp0] * diffL1[dd * (p + 1) + pp1] * L2[pp2] +
145  L0[pp0] * L1[pp1] * diffL2[dd * (p + 1) + pp2];
146  }
147  }
148  jj++;
149  }
150  }
151  }
152  }
153  if (jj != P)
154  SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P);
155  }
157 }
158 PetscErrorCode L2_VolumeShapeDiffMBTETinvJ(int base_p, int p,
159  double *volume_diffN, double *invJac,
160  double *volume_diffNinvJac,
161  int GDIM) {
163  int ii, gg;
164  for (ii = 0; ii != NBVOLUMETET_L2(p); ++ii) {
165  for (gg = 0; gg < GDIM; gg++) {
166  int shift1 = NBVOLUMETET_L2(base_p) * gg;
167  int shift2 = NBVOLUMETET_L2(p) * gg;
168  cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
169  &(volume_diffN)[3 * shift1 + 3 * ii], 1, 0.,
170  &(volume_diffNinvJac)[3 * shift2 + 3 * ii], 1);
171  }
172  }
174 }
175 
176 PetscErrorCode L2_ShapeFunctions_MBTRI(
177  int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM,
178  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
179  double *L, double *diffL,
180  const int dim)) {
181  return L2_Ainsworth_ShapeFunctions_MBTRI(p, N, diffN, L2N, diff_L2N, GDIM,
182  base_polynomials);
183 }
184 
185 PetscErrorCode L2_ShapeFunctions_MBTET(
186  int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM,
187  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
188  double *L, double *diffL,
189  const int dim)){
190  return L2_Ainsworth_ShapeFunctions_MBTET(p, N, diffN, L2N, diff_L2N, GDIM,
191  base_polynomials);
192 };
void cblas_dgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const double alpha, const double *A, const int lda, const double *X, const int incX, const double beta, double *Y, const int incY)
Definition: cblas_dgemv.c:11
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
static PetscErrorCode ierr
Definition: l2.c:27
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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:84
Loose implementation of some useful functions.
PetscErrorCode L2_VolumeShapeDiffMBTETinvJ(int base_p, int p, double *volume_diffN, double *invJac, double *volume_diffNinvJac, int GDIM)
Definition: l2.c:158
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
CHKERRQ(ierr)
useful compiler directives and definitions
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
Functions to approximate hierarchical spaces.
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:29
const int N
Definition: speed_test.cpp:3
PetscErrorCode L2_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))
Definition: l2.c:176
PetscErrorCode L2_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))
Definition: l2.c:185
field with C-1 continuity
Definition: definitions.h:174