add strong_boundaries

fluid/fluid_problem.c

@@ -506,12 +506,11 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
       double *local_vector = &all_local_vector[local_size*iel];
       double *local_matrix = &all_local_matrix[local_size*local_size*iel];
       const int *el = &mesh->elements[iel*N_N];
-      double dxidx[D][D], dphi[N_N][D];
+      double dxidx[D][D];
       mesh_dxidx(mesh,iel,dxidx);
       double n[D];
       double detbnd = get_normal(mesh, interface, n);
       const FEElement *element = mesh->element;
-      element->df(QP[0], dxidx, dphi);
       for (int iqp = 0; iqp < N_LQP; ++iqp) {
         double xiel[D];
         param_boundary_to_element(interface[1],LQP[iqp], xiel);
@@ -633,6 +632,9 @@ double fluid_problem_volume_flux(FluidProblem *problem, const char *tagname)
 
 double fluid_problem_a_integ_bnd(FluidProblem *problem, double dt)
 {
+
+
+
 #if 0
   const FEElement *element = problem->fields->element[0];
   const Mesh *mesh = problem->mesh;
@@ -1036,10 +1038,19 @@ void fluid_problem_assemble_system(FluidProblem *problem, const double *solution
   for (size_t i = 0; i < N_E*local_size*local_size; ++i)
     all_local_matrix[i] = 0;
 
-  // char *forced_row = malloc(sizeof(char)*n_fields*mesh->n_nodes);
   int ndof = fluid_problem_get_ndof(problem->fields, problem->mesh);
-  char *forced_row = malloc(sizeof(char)*ndof);
-  for (int i = 0; i < ndof; ++i)
+  fluid_problem_node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
+  compute_weak_boundary_conditions(problem, solution_old, dt, all_local_vector, all_local_matrix);
+  fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
+  if (problem->n_fluids == 2){
+    double *grad_a_cg = malloc(sizeof(double)*mesh->n_nodes*D);
+    fluid_problem_dg_to_cg_grad(problem, problem->concentration, problem->grad_a_cg);
+    fluid_problem_surface_tension(problem, solution_old, problem->grad_a_cg, all_local_vector);
+    free(grad_a_cg);
+  }
+
+  char *forced_row = malloc(sizeof(int)*ndof);
+  for(int i=0; i<ndof; ++i) 
     forced_row[i] = -1;
   for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
     const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
@@ -1052,40 +1063,49 @@ void fluid_problem_assemble_system(FluidProblem *problem, const double *solution
       printf("Boundary tag \"%s\" not found.\n", bnd->tag);
       exit(1);
     }
-    for (int i = 0; i < problem->boundaries[iphys].n_nodes; ++i){ // to change
-      forced_row[problem->boundaries[iphys].nodes[i]*n_fields+bnd->row] = bnd->field;
+    MeshBoundary *mbnd = &problem->boundaries[iphys];
+    const FEElement *field = problem->fields->element[bnd->field];
+    int ndof_closure = field->ndof_closure;
+    int c = 0;
+    for(int ifield = 0; ifield < problem->fields->n; ++ifield){
+        if (ifield == bnd->field)
+          break;
+        c += problem->fields->element[ifield]->nlocal;
+    }
+
+    for(int iinterface=0; iinterface < mbnd->n_interfaces; ++iinterface){
+      const int *interface = &mesh->interfaces[mbnd->interfaces[iinterface]*4];
+      const int iel = interface[0];
+      const int *el = &mesh->elements[iel*N_N];
+      const int *map = &(problem->fields->map[local_size*iel+c]);
+      for(int id=0; id < ndof_closure; ++id){
+        const int cl = field->closure[interface[1]*ndof_closure+id];
+        forced_row[map[cl]] = bnd->field;
+      }
     }
   }
-  fluid_problem_node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
-  compute_weak_boundary_conditions(problem, solution_old, dt, all_local_vector, all_local_matrix);
-  fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
-  if (problem->n_fluids == 2){
-    double *grad_a_cg = malloc(sizeof(double)*mesh->n_nodes*D);
-    fluid_problem_dg_to_cg_grad(problem, problem->concentration, problem->grad_a_cg);
-    fluid_problem_surface_tension(problem, solution_old, problem->grad_a_cg, all_local_vector);
-    free(grad_a_cg);
+  int c=0;
+  int c_phi[local_size], c_field[local_size];
+  for(int ifield=0; ifield<problem->fields->n; ++ifield){
+    for(int iphi=0; iphi < problem->fields->element[ifield]->nlocal; ++iphi){
+      c_field[c] = ifield;
+      c_phi[c] = iphi; 
+      c++;
+    }
+  }
+  for(int iel=0; iel < mesh->n_elements; ++iel){
+    double *local_vector = &all_local_vector[local_size*iel];
+    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
+    const int *map = &(problem->fields->map[local_size*iel]);
+    for(int ic = 0; ic < local_size; ++ic){
+        int forced_field = forced_row[map[ic]];
+        if (forced_field == -1) continue; 
+        local_vector[ic] = 0;
+        for(int jc=0 ; jc < local_size; ++jc){
+          local_matrix[ic*local_size+jc] = (c_phi[ic]==c_phi[jc] && c_field[jc]==forced_field)? 1. : 0.;
+        }
+    }
   }
-  // to do
-  // for (int iel=0; iel < mesh->n_elements; ++iel) {
-  //   double *local_vector = &all_local_vector[local_size*iel];
-  //   double *local_matrix = &all_local_matrix[local_size*local_size*iel];
-  //   for (int inode = 0; inode < element->nlocal; ++inode) {
-  //     for(int ifield = 0; ifield < n_fields; ++ifield) {
-  //       const int *el = &mesh->elements[iel*N_N];
-  //       // char forced_field = forced_row[el[inode]*n_fields+ifield];
-  //       char forced_field = -1;
-  //       if (forced_field == -1) continue;
-  //       int i = inode*n_fields + ifield;
-  //       for (int jnode = 0; jnode< element->nlocal; ++jnode){
-  //         for (int jfield= 0; jfield< n_fields; ++jfield){
-  //           local_matrix[(inode*n_fields+ifield)*local_size+(jnode*n_fields+jfield)] =
-  //             (inode==jnode && jfield == forced_field) ? 1. : 0.;
-  //         }
-  //       }
-  //       local_vector[inode+ifield*element->nlocal] = 0;
-  //     }
-  //   }
-  // }
   free(forced_row);
 }