fix some bugs

fluid/fluid_problem.c

@@ -328,7 +328,7 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
   free(h0);
 }
 
-static void f_pressure(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid, double *v){
+static void f_pressure(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid, double v){
   double p = s[P];
   double u[D], du[D][D], dp[D];
   for (int iD = 0; iD < D; ++iD) {
@@ -368,7 +368,7 @@ static void f_pressure(FluidProblem *problem,const double *n, double *f,const do
   }
 }
 
-static void f_velocity(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid,double *v){
+static void f_velocity(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid,double v){
   double p = s[P];
   double u[D], du[D][D], dp[D];
   for (int iD = 0; iD < D; ++iD) {
@@ -409,7 +409,7 @@ static void f_velocity(FluidProblem *problem,const double *n, double *f,const do
   }
 }
 
-static void f_wall(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds, const double c, const double *dc,const double dt,int eid,double *v)
+static void f_wall(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds, const double c, const double *dc,const double dt,int eid,double v)
 {
   double p = s[P];
   double u[D], du[D][D], dp[D];
@@ -607,6 +607,8 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
   int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}};
 #endif
   for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
+    for (int i =0; i < mesh->n_nodes; ++i)
+      bnd_node_local_id[i] = -1;
     if(mesh->phys_dimension[iphys] != D-1) continue;
     double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
     double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
@@ -633,8 +635,6 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
     }
 
     for (int ibnd = 0; ibnd < mesh->phys_n_boundaries[iphys]; ++ibnd) {
-      for (int i =0; i < mesh->n_nodes; ++i)
-        bnd_node_local_id[i] = -1;
       const int *bnd = &mesh->phys_boundaries[iphys][ibnd*2];
       const int iel = bnd[0];
       const int iebnd = bnd[1];
@@ -717,7 +717,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
           }
         }
         const double jw = LQW[iqp]*detbnd;
-        f(problem,n,f0,s,ds,c,dc,dt,iel,v);
+        f(problem,n,f0,s,ds,c,dc,dt,iel,vbnd);
         for (int ifield = 0; ifield < n_fields; ++ifield) {
           for (int iphi = 0; iphi < N_SF; ++iphi){
             local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
@@ -726,12 +726,12 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
         for (int jfield = 0; jfield < n_fields; ++jfield) {
           double store = s[jfield];
           s[jfield] += deps;
-          f(problem,n,g0,s,ds,c,dc,dt,iel,v);
+          f(problem,n,g0,s,ds,c,dc,dt,iel,vbnd);
           s[jfield] = store;
           for (int jd =0; jd < D; ++jd){
             store = ds[jfield*D+jd];
             ds[jfield*D+jd] += deps;
-            f(problem,n,h0+jd*n_fields,s,ds,c,dc,dt,iel,v);
+            f(problem,n,h0+jd*n_fields,s,ds,c,dc,dt,iel,vbnd);
             ds[jfield*D+jd] = store;
           }
           int N2 = n_fields*N_SF;

fluid/fluid_problem.h

@@ -43,7 +43,7 @@ typedef struct {
 } StrongBoundary;
 
 typedef struct FluidProblem FluidProblem;
-typedef void f_cb(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid, double *v);
+typedef void f_cb(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid, double v);
 typedef struct {
   char *tag;
   f_cb *cb;