keep track of the node tag and dofs for the dirichlet BCs

seamplex · Mar 1, 2025 · 1ff3dfd · 1ff3dfd
1 parent aaf7c42
commit 1ff3dfd
Show file tree

Hide file tree

Showing 13 changed files with 141 additions and 113 deletions.
diff --git a/configure.ac b/configure.ac
@@ -27,7 +27,7 @@ AC_PROG_INSTALL
 # default optimization flags
 AS_IF([test "x$CFLAGS" = "x-g -O2"], [CFLAGS="-O3 -flto=auto -no-pie"])
 # the CFLAGS are passed down to the linker as well
-# AS_IF([test "x$LDFLAGS" = "x"], [LDFLAGS="-flto=auto"])
+AS_IF([test "x$LDFLAGS" = "x"], [LDFLAGS="-flto=auto"])
 
 ######################
 # libraries
@@ -391,11 +391,10 @@ AC_DEFINE_UNQUOTED([FEENOX_COMPILER_LDFLAGS], ["${LDFLAGS}"], [LDFLAGS])
 # create some links to pass make distcheck
 m4_include(auto_links.m4)
 
-
-#AC_CONFIG_FILES([Makefile src/Makefile src/benchmark/Makefile doc/Makefile]) 
-#AC_CONFIG_FILES([Makefile src/Makefile doc/Makefile])
-#AC_CONFIG_FILES([Makefile src/Makefile]) 
-AC_CONFIG_FILES([Makefile src/Makefile utils/fee2ccx/Makefile]) 
+AC_CONFIG_FILES([Makefile src/Makefile])
+# AC_CONFIG_FILES([src/benchmark/Makefile]) 
+# AC_CONFIG_FILES([doc/Makefile])
+AS_IF([test -d utils/fee2ccx],[AC_CONFIG_FILES([utils/fee2ccx/Makefile])]) 
 
 
 AS_BOX([Summary of dependencies])

diff --git a/src/feenox.h b/src/feenox.h
@@ -1923,6 +1923,8 @@ struct feenox_t {
 
 
     // internal storage of evaluated dirichlet conditions
+    size_t          *dirichlet_nodes;
+    unsigned int    *dirichlet_dofs;
     PetscInt        *dirichlet_indexes;
     PetscScalar     *dirichlet_values;
     PetscScalar     *dirichlet_derivatives;
@@ -2401,7 +2403,7 @@ extern PetscErrorCode feenox_ts_jacobian(TS ts, PetscReal t, Vec T, Vec T_dot, P
 extern int feenox_problem_solve_slepc_eigen(void);
 
 // dirichlet.c
-extern int feenox_problem_dirichlet_add(size_t index, double value);
+extern int feenox_problem_dirichlet_add(size_t j_global, unsigned int g, double value);
 extern int feenox_problem_multifreedom_add(size_t index, double *coefficients);
 extern int feenox_problem_mimicked_add(size_t index_guide, size_t index_follower, double coefficient_follower);
 extern int feenox_problem_dirichlet_eval(void);

diff --git a/src/math/assignment.c b/src/math/assignment.c
@@ -178,7 +178,7 @@ int feenox_add_assignment(const char *left_hand, const char *right_hand) {
         feenox_push_error_message("right hand side for vector initialization has to be inside brackets '(' and ')'");
         return FEENOX_ERROR;
       }
-      *rhs_starting_with_bracket++;
+      rhs_starting_with_bracket++;
     }
     size_t n_chars_count = 0;
     int n_expressions = 0;

diff --git a/src/mesh/calculix.c b/src/mesh/calculix.c
@@ -130,7 +130,7 @@ int frdfromgmsh_types[] = {
  ELEMENT_TYPE_PRISM6,        // 2
  ELEMENT_TYPE_TETRAHEDRON4,  // 3
  ELEMENT_TYPE_HEXAHEDRON20,  // 4
- ELEMENT_TYPE_UNDEFINED,     // 5
+ ELEMENT_TYPE_PRISM15,       // 6
  ELEMENT_TYPE_TETRAHEDRON10, // 6
  ELEMENT_TYPE_TRIANGLE3,     // 7
  ELEMENT_TYPE_TRIANGLE6,     // 8

diff --git a/src/pdes/dirichlet.c b/src/pdes/dirichlet.c
@@ -21,13 +21,16 @@
  */
 #include "feenox.h"
 
-int feenox_problem_dirichlet_add(size_t index, double value) {
+int feenox_problem_dirichlet_add(size_t j_global, unsigned int g, double value) {
 
+  feenox.pde.dirichlet_nodes[feenox.pde.dirichlet_k] = feenox.pde.mesh->node[j_global].tag;
+  feenox.pde.dirichlet_dofs[feenox.pde.dirichlet_k] = g;
+
 #ifdef HAVE_PETSC
-  feenox.pde.dirichlet_indexes[feenox.pde.dirichlet_k] = index;
+  feenox.pde.dirichlet_indexes[feenox.pde.dirichlet_k] = feenox.pde.mesh->node[j_global].index_dof[g];
   feenox.pde.dirichlet_values[feenox.pde.dirichlet_k] = value;
-  feenox.pde.dirichlet_k++;
 #endif
+  feenox.pde.dirichlet_k++;
 
   return FEENOX_OK;
 }
@@ -69,6 +72,8 @@ int feenox_problem_dirichlet_eval(void) {
     // TODO: allow the user to provide a factor at runtime
     n_bcs = 1.2*feenox.pde.dofs * (feenox.pde.last_node - feenox.pde.first_node);    
 
+    feenox_check_alloc(feenox.pde.dirichlet_nodes = calloc(n_bcs, sizeof(size_t)));
+    feenox_check_alloc(feenox.pde.dirichlet_dofs = calloc(n_bcs, sizeof(int)));
     feenox_check_alloc(feenox.pde.dirichlet_indexes = calloc(n_bcs, sizeof(PetscInt)));
     feenox_check_alloc(feenox.pde.dirichlet_values = calloc(n_bcs, sizeof(PetscScalar)));
     feenox_check_alloc(feenox.pde.dirichlet_derivatives = calloc(n_bcs, sizeof(PetscScalar)));
@@ -119,6 +124,8 @@ int feenox_problem_dirichlet_eval(void) {
     feenox.pde.dirichlet_rows = feenox.pde.dirichlet_k;
 
     // if k == 0 this like freeing
+    feenox_check_alloc(feenox.pde.dirichlet_nodes = realloc(feenox.pde.dirichlet_nodes, feenox.pde.dirichlet_rows * sizeof(size_t)));
+    feenox_check_alloc(feenox.pde.dirichlet_dofs = realloc(feenox.pde.dirichlet_dofs, feenox.pde.dirichlet_rows * sizeof(int)));
     feenox_check_alloc(feenox.pde.dirichlet_indexes = realloc(feenox.pde.dirichlet_indexes, feenox.pde.dirichlet_rows * sizeof(PetscInt)));
     feenox_check_alloc(feenox.pde.dirichlet_values = realloc(feenox.pde.dirichlet_values, feenox.pde.dirichlet_rows * sizeof(PetscScalar)));
     feenox_check_alloc(feenox.pde.dirichlet_derivatives = realloc(feenox.pde.dirichlet_derivatives, feenox.pde.dirichlet_rows * sizeof(PetscScalar)));
@@ -140,8 +147,7 @@ int feenox_problem_dirichlet_eval(void) {
 // it takes K and b and writes K_bc and b_bc
 int feenox_problem_dirichlet_set_K(void) {
 
-  if (feenox.pde.dirichlet_scale == 0)
-  {
+  if (feenox.pde.dirichlet_scale == 0) {
     feenox_problem_dirichlet_compute_scale();
   }
 

diff --git a/src/pdes/laplace/bc.c b/src/pdes/laplace/bc.c
@@ -61,7 +61,7 @@ int feenox_problem_bc_set_laplace_phi(bc_data_t *this, element_t *e, size_t j_gl
 
 #ifdef HAVE_PETSC
 
-  feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[0], feenox_expression_eval(&this->expr)));
+  feenox_call(feenox_problem_dirichlet_add(j_global, 0, feenox_expression_eval(&this->expr)));
   // TODO: only in transient
 //  feenox.pde.dirichlet_derivatives[*k] = feenox_expression_derivative_wrt_variable(&bc_data->expr, feenox_special_var(t), feenox_special_var_value(t));
 

diff --git a/src/pdes/mechanical/bc.c b/src/pdes/mechanical/bc.c
@@ -156,16 +156,16 @@ int feenox_problem_bc_set_mechanical_displacement(bc_data_t *this, element_t *e,
   if (this->dof != -1) {
 
     // only one dof
-    feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[this->dof], feenox_expression_eval(&this->expr)));
+    feenox_call(feenox_problem_dirichlet_add(j_global, this->dof, feenox_expression_eval(&this->expr)));
 
   } else {
 
     // -1 means all dofs (and the only possibility is to have all them equal to zero)
-    feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[0], 0));
+    feenox_call(feenox_problem_dirichlet_add(j_global, 0, 0));
     if (feenox.pde.dofs > 1) {
-      feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[1], 0));
+      feenox_call(feenox_problem_dirichlet_add(j_global, 1, 0));
       if (feenox.pde.dofs > 2) {
-        feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[2], 0));
+        feenox_call(feenox_problem_dirichlet_add(j_global, 2, 0));
       }
     }
 
@@ -213,7 +213,7 @@ int feenox_problem_bc_set_mechanical_symmetry(bc_data_t *this, element_t *e, siz
   // then we set a traditional dirichlet bc (i.e. u=0 or v=0 or w=0)
   // otherwise we need a generic multifreedom
   if (coordinate_direction != -1) {
-    feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[coordinate_direction], 0));
+    feenox_call(feenox_problem_dirichlet_add(j_global, coordinate_direction, 0));
   } else {
     feenox_call(feenox_problem_multifreedom_add(j_global, normal));
   }  

diff --git a/src/pdes/modal/bc.c b/src/pdes/modal/bc.c
@@ -117,16 +117,16 @@ int feenox_problem_bc_set_modal_displacement(bc_data_t *this, element_t *e, size
   if (this->dof != -1) {
 
     // only one dof
-    feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[this->dof], feenox_expression_eval(&this->expr)));
+    feenox_call(feenox_problem_dirichlet_add(j_global, this->dof, feenox_expression_eval(&this->expr)));
 
   } else {
 
     // -1 means all dofs (and the only possibility is to have all them equal to zero)
-    feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[0], 0));
+    feenox_call(feenox_problem_dirichlet_add(j_global, 0, 0));
     if (feenox.pde.dofs > 1) {
-      feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[1], 0));
+      feenox_call(feenox_problem_dirichlet_add(j_global, 1, 0));
       if (feenox.pde.dofs > 2) {
-        feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[2], 0));
+        feenox_call(feenox_problem_dirichlet_add(j_global, 2, 0));
       }
     }
 
@@ -177,7 +177,7 @@ int feenox_problem_bc_set_modal_symmetry(bc_data_t *this, element_t *e, size_t j
   // then we set a traditional dirichlet bc (i.e. u=0 or v=0 or w=0)
   // otherwise we need a generic multifreedom
   if (coordinate_direction != -1) {
-    feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[coordinate_direction], 0));
+    feenox_call(feenox_problem_dirichlet_add(j_global, coordinate_direction, 0));
   } else {
     feenox_call(feenox_problem_multifreedom_add(j_global, normal));
   }  

diff --git a/src/pdes/neutron_diffusion/bc.c b/src/pdes/neutron_diffusion/bc.c
@@ -106,7 +106,7 @@ int feenox_problem_bc_set_neutron_diffusion_null(bc_data_t *this, element_t *e,
 
 #ifdef HAVE_PETSC
   for (unsigned int g = 0; g < feenox.pde.dofs; g++) {
-    feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[g], 0));
+    feenox_call(feenox_problem_dirichlet_add(j_global, g, 0));
   }
 #endif
 
@@ -118,7 +118,7 @@ int feenox_problem_bc_set_neutron_diffusion_null(bc_data_t *this, element_t *e,
 int feenox_problem_bc_set_neutron_diffusion_flux(bc_data_t *this, element_t *e, size_t j_global) {
 
 #ifdef HAVE_PETSC
-  feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[this->dof], feenox_expression_eval(&this->expr)));
+  feenox_call(feenox_problem_dirichlet_add(j_global, this->dof, feenox_expression_eval(&this->expr)));
 #endif
   return FEENOX_OK;
 }

diff --git a/src/pdes/neutron_sn/bc.c b/src/pdes/neutron_sn/bc.c
@@ -70,7 +70,7 @@ int feenox_problem_bc_set_neutron_sn_vacuum(bc_data_t *this, element_t *e, size_
     if (feenox_mesh_dot(neutron_sn.Omega[m], outward_normal) < 0) {
       // if the direction is inward set it to zero
       for (unsigned int g = 0; g < neutron_sn.groups; g++) {
-        feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[sn_dof_index(m,g)], 0));
+        feenox_call(feenox_problem_dirichlet_add(j_global, sn_dof_index(m,g), 0));
       }
     }
   }

diff --git a/src/pdes/thermal/bc.c b/src/pdes/thermal/bc.c
@@ -88,7 +88,7 @@ int feenox_problem_bc_set_thermal_temperature(bc_data_t *this, element_t *e, siz
 
 #ifdef HAVE_PETSC
 
-  feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[0], feenox_expression_eval(&this->expr)));  
+  feenox_call(feenox_problem_dirichlet_add(j_global, 0, feenox_expression_eval(&this->expr)));  
   // TODO: only in transient
 //  feenox.pde.dirichlet_derivatives[*k] = feenox_expression_derivative_wrt_variable(&bc_data->expr, feenox_special_var(t), feenox_special_var_value(t));
 

diff --git a/utils/fee2ccx/configure.ac b/utils/fee2ccx/configure.ac