wip

fluid/fluid_problem_fem.c

@@ -63,8 +63,22 @@ static void fe_df_triangle_p2(const double *xi, const double dxidx[][D], double
   }
 }
 
+static void fe_f_tetrahedron_p1(const double *xi, double *f) {
+  f[0] = 1-xi[0]-xi[1]-xi[2];
+  f[1] = xi[0];
+  f[2] = xi[1];
+  f[3] = xi[2];
+}
 
+static void fe_df_tetrahedron_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]-dxidx[2][d];
+    df[1][d] = dxidx[0][d];
+    df[2][d] = dxidx[1][d];
+    df[3][d] = dxidx[2][d];
+  }
 
+}
 FEElement *fe_element_new(const char *etype) {
   FEElement *fe = malloc(sizeof(FEElement));
   if (strcmp(etype,"triangle_p1") == 0) {
@@ -112,6 +126,16 @@ FEElement *fe_element_new(const char *etype) {
     fe->f = fe_f_line_p2;
     fe->df = fe_df_line_p2;
   }
+  else if (strcmp(etype,"tetrahedron_p1") == 0) {
+    fe->nlocal = 4;
+    fe->n[0] = 1;
+    fe->n[1] = 0;
+    fe->n[2] = 0;
+    fe->n[3] = 0;
+    fe->f = fe_f_tetrahedron_p1;
+    fe->df = fe_df_tetrahedron_p1;
+  }
+
   else {
     printf("Unknown element type '%s'.\n", etype);
     exit(0);