wip

fluid/fluid_problem.c

@@ -109,16 +109,16 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
     double *local_matrix = all_local_matrix ? &all_local_matrix[fields->local_size*fields->local_size*iel] : NULL;
     double c, cold, dc[D];
     double s[n_fields], ds[n_fields][D];
-    fe_fields_eval_grad(fields, mesh, iel, sf, dsf, problem->solution, s, ds);
+    fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, problem->solution, s, ds);
     double sold[n_fields];
-    fe_fields_eval(fields, mesh, iel, sf,solution_old, sold);
-    fe_fields_eval_grad(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
-    fe_fields_eval(problem->field_porosity, mesh, iel, sfporosity, problem->oldporosity, &cold);
+    fe_fields_eval_sf(fields, mesh, iel, sf,solution_old, sold);
+    fe_fields_eval_grad_sf(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
+    fe_fields_eval_sf(problem->field_porosity, mesh, iel, sfporosity, problem->oldporosity, &cold);
     double a = 0;
     if (problem->n_fluids == 2) {
       double sfconcentration[problem->field_concentration->local_size];
       fe_fields_f(problem->field_concentration, xi, sfconcentration);
-      fe_fields_eval(problem->field_concentration, mesh, iel, problem->concentration, sfconcentration, &a);
+      fe_fields_eval_sf(problem->field_concentration, mesh, iel, problem->concentration, sfconcentration, &a);
     }
     double f[D],dfdu,dfddp;
     particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
@@ -513,12 +513,12 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
         double sf[fields->local_size], dsf[fields->local_size][D];
         fe_fields_f(fields, xiel, sf);
         fe_fields_df(fields, QP[iqp], dxidx, dsf);
-        fe_fields_eval_grad(fields, mesh, iel, sf, dsf, problem->solution, s, ds);
-        fe_fields_eval_grad(fields, mesh, iel, sf, dsf, solution_old, sold, dsold);
+        fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, problem->solution, s, ds);
+        fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, solution_old, sold, dsold);
         double dc[D], c, sfporosity[problem->field_porosity->local_size], dsfporosity[problem->field_porosity->local_size][D] ;
         fe_fields_f(problem->field_porosity, QP[iqp], sfporosity);
         fe_fields_df(problem->field_porosity, QP[iqp], dxidx, dsfporosity);
-        fe_fields_eval_grad(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
+        fe_fields_eval_grad_sf(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
         double f0[fields->n], f00[fields->n*fields->n],f01[fields->n*fields->n][D];
         for (int i = 0; i < fields->n; ++i) {
           f0[i] = 0;
@@ -535,7 +535,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
         if (problem->n_fluids==2) {
           double sfconcentration[problem->field_concentration->local_size];
           fe_fields_f(problem->field_concentration, QP[iqp], sfconcentration);
-          fe_fields_eval(problem->field_concentration, mesh, iel, sfconcentration, problem->concentration, &a);
+          fe_fields_eval_sf(problem->field_concentration, mesh, iel, sfconcentration, problem->concentration, &a);
         }
         const FEElement *mesh_element = problem->mesh->element;
         double meshf[mesh_element->nlocal];
@@ -731,19 +731,19 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
       double sf[fields->local_size], dsf[fields->local_size][D];
       fe_fields_f(fields, QP[iqp], sf);
       fe_fields_df(fields, QP[iqp], dxidx, dsf);
-      fe_fields_eval_grad(fields, mesh, iel, sf, dsf, solution, s, ds);
-      fe_fields_eval_grad(fields, mesh, iel, sf, dsf, solution_old, sold, dsold);
+      fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, solution, s, ds);
+      fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, solution_old, sold, dsold);
       double sfporosity[fields->local_size], dsfporosity[fields->local_size][D];
       fe_fields_f(problem->field_porosity, QP[iqp], sfporosity);
       fe_fields_df(problem->field_porosity, QP[iqp], dxidx, dsfporosity);
       double c, dc[D], cold;
-      fe_fields_eval_grad(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
-      fe_fields_eval(problem->field_porosity, mesh, iel, sfporosity, problem->porosity, &cold);
+      fe_fields_eval_grad_sf(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
+      fe_fields_eval_sf(problem->field_porosity, mesh, iel, sfporosity, problem->porosity, &cold);
       double a = 0;
       if (problem->n_fluids==2) {
         double sfconcentration[problem->field_concentration->local_size];
         fe_fields_f(problem->field_concentration, QP[iqp], sfconcentration);
-        fe_fields_eval(problem->field_concentration, mesh, iel, sfconcentration, problem->concentration, &a);
+        fe_fields_eval_sf(problem->field_concentration, mesh, iel, sfconcentration, problem->concentration, &a);
       }
       double bf[D];
       const FEElement *mesh_element = problem->mesh->element;

fluid/fluid_problem.h

@@ -101,7 +101,7 @@ struct FluidProblem {
   int advection;
   int temporal;
   int drag_in_stab;
-  FEFields *fields, *field_porosity, *field_concentration;
+  FEFields *fields, *field_porosity, *field_concentration, *fields_bulk_force;
 };
 
 // complete force on the particle (including gravity)

fluid/fluid_problem_fem.h

@@ -15,18 +15,20 @@ struct FEElementStruct{
 FEElement *fe_element_new(const char *etype);
 void fe_element_free(FEElement *element);
 
-typedef struct {
+struct FEFieldsStruct{
   int local_size;
   int n;
   int nd[DIMENSION+1];
   FEElement **element;
-}FEFields;
+};
 
 FEFields *fe_fields_new(const char **etypes);
 void fe_fields_f(const FEFields *fields, const double *xi, double *f);
 void fe_fields_df(const FEFields *fields, const double *xi, const double dxdidx[D][D], double df[][D]);
-void fe_fields_eval_grad(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double dsf[][DIMENSION], const double *v, double s[], double ds[][DIMENSION]);
-void fe_fields_eval(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double *v, double s[]);
+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]);
+void fe_fields_eval_sf(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double *v, double s[]);
+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]);
+void fe_fields_eval(const FEFields *fields, const double *xi, const Mesh *mesh, int iel, const double *v, double s[]);
 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);
 static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *f00, const double f01[][D], const double f10[][D], const double f11[][D][D], const double *sf, const double dsf[][D], double jw, double *local_matrix) {
   int jc = 0;

fluid/mesh.c

@@ -154,6 +154,7 @@ static void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *
 
 void mesh_free(Mesh *m)
 {
+  fe_fields_free(m->coord_fields);
   free(m->elements);
   free(m->x);
   for (int i = 0; i < m->n_boundaries; ++i) {
@@ -173,6 +174,13 @@ Mesh *mesh_new_from_elements(int n_nodes, double *x, int n_elements, int *elemen
 #else
   m->element = fe_element_new("tetrahedron_p1");
 #endif
+#if DIMENSION == 3
+  const char *types[] = {"tetrahedron_p1","tetrahedron_p1", "tetrahedron_p1", NULL};
+#else
+  const char *types[] = {"triangle_p1","triangle_p1",NULL};
+#endif
+  m->coord_fields = fe_fields_new(types);
+
   m->n_elements = n_elements;
   m->elements = malloc(n_elements*sizeof(int)*(DIMENSION+1));
   memcpy(m->elements,elements,(DIMENSION+1)*sizeof(int)*n_elements);

fluid/mesh.h

@@ -27,6 +27,7 @@
 #include <stdio.h>
 #include <stdint.h>
 typedef struct FEElementStruct FEElement;
+typedef struct FEFieldsStruct FEFields;
 
 typedef struct {
   int n_nodes;
@@ -42,6 +43,7 @@ typedef struct {
   int n_periodic;
   int *periodic_mapping; // n_nodes
   FEElement *element;
+  FEFields *coord_fields;
 } Mesh;
 
 void mesh_free(Mesh *m);