missing file

fluid/fluid_problem_fem.c

+#include "fluid_problem_fem.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+static void fe_f_line_p1(const double *xi, double *f) {
+  f[0] = 1-xi[0];
+  f[1] = xi[0];
+}
+static void fe_df_line_p1(const double *xi, const double dxidx[][D], double df[][D]) {
+  for (int d = 0; d < D; ++d) {
+    df[0][d] = -dxidx[0][d];
+    df[1][d] = dxidx[0][d];
+  }
+}
+
+static void fe_f_triangle_p1(const double *xi, double *f) {
+  f[0] = 1-xi[0]-xi[1];
+  f[1] = xi[0];
+  f[2] = xi[1];
+}
+
+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];
+    df[2][d] = dxidx[1][d];
+  }
+}
+
+
+FEElement *fe_element_new(const char *etype) {
+  FEElement *fe = malloc(sizeof(FEElement));
+  if (strcmp(etype,"triangle_p1") == 0) {
+    fe->nlocal = 3;
+    fe->n[0] = 1;
+    fe->n[1] = 0;
+    fe->n[2] = 0;
+    fe->n[3] = 0;
+    fe->f = fe_f_triangle_p1;
+    fe->df = fe_df_triangle_p1;
+  }
+  else if (strcmp(etype,"triangle_p1dg") == 0) {
+    fe->nlocal = 3;
+    fe->n[0] = 0;
+    fe->n[1] = 0;
+    fe->n[2] = 3;
+    fe->n[3] = 0;
+    fe->f = fe_f_triangle_p1;
+    fe->df = fe_df_triangle_p1;
+  }
+  else if (strcmp(etype,"line_p1") == 0) {
+    fe->nlocal = 2;
+    fe->n[0] = 1;
+    fe->n[1] = 0;
+    fe->n[2] = 0;
+    fe->n[3] = 0;
+    fe->f = fe_f_line_p1;
+    fe->df = fe_df_line_p1;
+  }
+  else {
+    printf("Unknown element type '%s'.\n", etype);
+    exit(0);
+  }
+  return fe;
+}
+
+void fe_element_free(FEElement *fe) {
+  free(fe);
+}
+
+
+FEFields *fe_fields_new(const char **etypes) {
+  FEFields *fields = malloc(sizeof(FEFields));
+  fields->n = 0;
+  fields->local_size = 0;
+  for (const char **etype = etypes; *etype != NULL; etype++) {
+    fields->n++;
+  }
+  fields->element = malloc(sizeof(FEElement*)*fields->n);
+  int n = 0;
+  for (int i = 0; i < fields->n; ++i) {
+    fields->element[i] = fe_element_new(etypes[i]);
+    fields->local_size += fields->element[i]->nlocal;
+  }
+  for (int d = 0; d < DIMENSION+1; ++d) {
+    fields->nd[d] = 0;
+    for (int i = 0; i < fields->n; ++i) {
+      fields->nd[d] += fields->element[i]->n[d];
+    }
+  }
+  return fields;
+}
+
+void fe_fields_free(FEFields *fields) {
+  for (int i = 0; i < fields->n; ++i) {
+    fe_element_free(fields->element[i]);
+  }
+  free(fields->element);
+  free(fields);
+}
+
+void fe_fields_f(const FEFields *fields, const double *xi, double *f) {
+  int c = 0;
+  for (int i = 0; i < fields->n; ++i) {
+    fields->element[i]->f(xi, &f[c]);
+    c += fields->element[i]->nlocal;
+  }
+}
+
+void fe_fields_df(const FEFields *fields, const double *xi, const double dxidx[D][D], double df[][D]) {
+  int c = 0;
+  for (int i = 0; i < fields->n; ++i) {
+    fields->element[i]->df(xi, dxidx, &df[c]);
+    c += fields->element[i]->nlocal;
+  }
+}
+
+void fe_fields_local_map(const FEFields *fields, const Mesh *mesh, int iel, int map[]) {
+  int c[fields->n];
+  c[0] = 0;
+  for (int i = 1; i < fields->n; ++i) {
+    c[i] = c[i-1]+fields->element[i-1]->nlocal;
+  }
+  int l = 0;
+  for (int i = 0; i < fields->n; ++i) {
+    for (int j = 0; j < fields->element[i]->n[0]; ++j){
+      for (int k = 0; k < DIMENSION+1; ++k) {
+        map[c[i]++] = mesh->elements[iel*(DIMENSION+1)+k]*fields->nd[0]+l;
+      }
+      l++;
+    }
+  }
+  int shift = l*mesh->n_nodes;
+  // todo add edges (and faces in 3D)
+  for (int i = 0; i < fields->n; ++i) {
+    for (int j = 0; j < fields->element[i]->n[DIMENSION]; ++j){
+      map[c[i]++] = shift+iel*fields->nd[DIMENSION]+l;
+      l++;
+    }
+  }
+  if (fields->n == 3) {
+    int *el = mesh->elements + 3*iel;
+    int map2[9] = {
+      el[0]*3+0, el[1]*3+0, el[2]*3+0,
+      el[0]*3+1, el[1]*3+1, el[2]*3+1,
+      el[0]*3+2, el[1]*3+2, el[2]*3+2,
+    };
+  }
+}
+
+void fe_fields_eval_sf(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double *v, double s[]) {
+  int map[fields->local_size];
+  fe_fields_local_map(fields, mesh, iel, map);
+  int c = 0;
+  for (int i = 0; i < fields->n; ++i) {
+    s[i] = 0;
+    for (int j = 0; j < fields->element[i]->nlocal; ++j) {
+      s[i] += v[map[c]]*sf[c];
+      c++;
+    }
+  }
+}
+
+void fe_fields_eval_grad_sf(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double dsf[][DIMENSION], const double *v, double s[], double ds[][DIMENSION]) {
+  int map[fields->local_size];
+  fe_fields_local_map(fields, mesh, iel, map);
+  int c = 0;
+  for (int i = 0; i < fields->n; ++i) {
+    s[i] = 0;
+    for (int d = 0; d < DIMENSION; ++d) {
+      ds[i][d] = 0;
+    }
+    for (int j = 0; j < fields->element[i]->nlocal; ++j) {
+      s[i] += v[map[c]]*sf[c];
+      for (int d = 0; d < DIMENSION; ++d){
+        ds[i][d] += v[map[c]]*dsf[c][d];
+      }
+      c++;
+    }
+  }
+}
+
+void fe_fields_add_to_local_vector(const FEFields *fields, const double *f0, const double f1[][D], const double *sf, const double dsf[][D], double jw, double *local_vector) {
+  int cc = 0;
+  for (int ifield = 0; ifield < fields->n; ++ifield) {
+    for (int iphi = 0; iphi < fields->element[ifield]->nlocal; ++iphi){
+      double v = 0;
+      if (f0) {
+        v += sf[cc]*f0[ifield]*jw;
+      }
+      if (f1){
+        for (int id = 0; id < D; ++id) {
+          v += dsf[cc][id]*f1[ifield][id]*jw;
+        }
+      }
+      local_vector[cc] += v;
+      cc++;
+    }
+  }
+}
+
+void fe_fields_eval_grad(const FEFields *fields, const double *xi, const double dxidx[D][D], const Mesh *mesh, int iel, const double *v, double s[], double ds[][DIMENSION]) {
+  double sf[fields->local_size], dsf[fields->local_size][D];
+  fe_fields_f(fields, xi, sf);
+  fe_fields_df(fields, xi, dxidx, dsf);
+  fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, v, s, ds);
+}
+
+void fe_fields_eval(const FEFields *fields, const double *xi, const Mesh *mesh, int iel, const double *v, double s[]) {
+  double sf[fields->local_size];
+  fe_fields_f(fields, xi, sf);
+  fe_fields_eval_sf(fields, mesh, iel, sf, v, s);
+}