triangle_p2

fluid/fluid_problem_fem.c

@@ -20,7 +20,7 @@ static void fe_f_triangle_p1(const double *xi, double *f) {
   f[2] = xi[1];
 }
 
-void fe_df_triangle_p1(const double *xi, const double dxidx[][D], double df[][D]) {
+static void fe_df_triangle_p1(const double *xi, const double dxidx[][D], double df[][D]) {
   for (int d = 0; d < D; ++d) {
     df[0][d] = -dxidx[0][d]-dxidx[1][d];
     df[1][d] = dxidx[0][d];
@@ -28,6 +28,27 @@ void fe_df_triangle_p1(const double *xi, const double dxidx[][D], double df[][D]
   }
 }
 
+static void fe_f_triangle_p2(const double *xi, double *f) {
+  f[0] = 1-3*(xi[0]+xi[1])+2*(xi[0]+xi[1])*(xi[0]+xi[1]);
+  f[1] = xi[0]*(2*xi[0]-1);
+  f[2] = xi[1]*(2*xi[1]-1);
+  f[3] = 4*xi[0]*(1-xi[0]-xi[1]);
+  f[4] = 4*xi[0]*xi[1];
+  f[5] = 4*xi[1]*(1-xi[0]-xi[1]);
+}
+
+static void fe_df_triangle_p2(const double *xi, const double dxidx[][D], double df[][D]) {
+  for (int d = 0; d < D; ++d) {
+    df[0][d] = (-3+4*(xi[0]+xi[1]))*(dxidx[0][d] + dxidx[1][d]);
+    df[1][d] = (4*xi[0]-1)*dxidx[0][d];
+    df[2][d] = (4*xi[1]-1)*dxidx[1][d];
+    df[3][d] = 4*((1-2*xi[0]-xi[1])*dxidx[0][d]-xi[0]*dxidx[1][d]);
+    df[4][d] = 4*(xi[1]*dxidx[0][d] + xi[0]*dxidx[1][d]);
+    df[5][d] = 4*(-xi[1]*dxidx[0][d] + (1-xi[0]-2*xi[1])*dxidx[1][d]);
+  }
+}
+
+
 
 FEElement *fe_element_new(const char *etype) {
   FEElement *fe = malloc(sizeof(FEElement));
@@ -49,6 +70,15 @@ FEElement *fe_element_new(const char *etype) {
     fe->f = fe_f_triangle_p1;
     fe->df = fe_df_triangle_p1;
   }
+  if (strcmp(etype,"triangle_p2") == 0) {
+    fe->nlocal = 6;
+    fe->n[0] = 1;
+    fe->n[1] = 1;
+    fe->n[2] = 0;
+    fe->n[3] = 0;
+    fe->f = fe_f_triangle_p2;
+    fe->df = fe_df_triangle_p2;
+  }
   else if (strcmp(etype,"line_p1") == 0) {
     fe->nlocal = 2;
     fe->n[0] = 1;
@@ -139,6 +169,7 @@ void fe_fields_local_map(const FEFields *fields, const Mesh *mesh, int iel, int
       l++;
     }
   }
+  // unused
   if (fields->n == 3) {
     int *el = mesh->elements + 3*iel;
     int map2[9] = {