Investigation of segfault with O3 when computing outward normal

Makefile.am

@@ -19,7 +19,6 @@ dist_doc_DATA = AUTHORS \
 
 dist_man_MANS = doc/feenox.1
 
-
 TESTS = \
         tests/abort.sh \
         tests/algebraic_expr.sh \

configure.ac

@@ -25,7 +25,7 @@ AC_PROG_INSTALL
 
 ######################
 # default optimization flags
-AS_IF([test "x$CFLAGS" = "x-g -O2"], [CFLAGS="-O2 -flto=auto -no-pie"])
+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"])
 

src/feenox.h

@@ -2314,7 +2314,6 @@ extern int feenox_mesh_elements_info(void);
 extern int feenox_mesh_elements_info_geo(element_type_t *element_type);
 
 // geom.c
-// TODO: rename mesh_* -> feenox_mesh_*
 extern void feenox_mesh_subtract(const double *a, const double *b, double *c);
 extern void feenox_mesh_cross(const double *a, const double *b, double *c);
 extern void feenox_mesh_normalized_cross(const double *a, const double *b, double *c);
@@ -2327,7 +2326,7 @@ extern double feenox_mesh_subtract_module(const double *b, const double *a);
 extern double feenox_mesh_subtract_squared_module(const  double *b, const  double *a);
 extern double feenox_mesh_subtract_squared_module2d(const  double *b, const  double *a);
 
-extern int feenox_mesh_compute_outward_normal(element_t *element, double n[3]);
+extern int feenox_mesh_compute_outward_normal(element_t *element, double n[3]) __attribute__((noinline));
 
 // element.c
 extern int feenox_mesh_compute_normal_2d(element_t *e);
@@ -2348,8 +2347,8 @@ extern int feenox_mesh_write_data_vtk(mesh_write_t *, mesh_write_dist_t *dist);
 
 
 // neighbors.c
-extern element_t *feenox_mesh_find_element_volumetric_neighbor(element_t *);
-extern int feenox_mesh_count_element_volumetric_neighbors(element_t *this);
+extern element_t *feenox_mesh_find_element_volumetric_neighbor(element_t *e)  __attribute__((noinline));
+extern int feenox_mesh_count_element_volumetric_neighbors(element_t *e)  __attribute__((noinline));
 
 
 // init.c

src/mesh/geometry.c

@@ -92,33 +92,33 @@ double feenox_mesh_subtract_squared_module2d(const  double *b, const  double *a)
   return (b[0]-a[0])*(b[0]-a[0]) + (b[1]-a[1])*(b[1]-a[1]);
 }
 
-
 // TODO: make a faster one assuming the elements are already oriented
 int feenox_mesh_compute_outward_normal(element_t *element, double n[3]) {
 
-  double local_n[3] = {n[0], n[1], n[2]};
-
+  // double local_n[3] = {n[0], n[1], n[2]};
+  // element_t *local_element = calloc(1, sizeof(element_t));
+  // memcpy(local_element, element, sizeof(element_t));
   // TODO: a method linked to the element type
   if (element->type->dim == 0) {
-    local_n[0] = 1;
-    local_n[1] = 0;
-    local_n[2] = 0;
-    
+    n[0] = 1;
+    n[1] = 0;
+    n[2] = 0;
+
   } else if (element->type->dim == 1) {
 
     // WATCH out! this does not work with lines which do not lie on the xy plane!
     double module = feenox_mesh_subtract_module(element->node[1]->x, element->node[0]->x);
-    local_n[0] = -(element->node[1]->x[1] - element->node[0]->x[1])/module;
-    local_n[1] = +(element->node[1]->x[0] - element->node[0]->x[0])/module;
-    local_n[2] = 0;
+    n[0] = -(element->node[1]->x[1] - element->node[0]->x[1])/module;
+    n[1] = +(element->node[1]->x[0] - element->node[0]->x[0])/module;
+    n[2] = 0;
 
   } else if (element->type->dim == 2) {
 
     double a[3] = {0,0,0};
     double b[3] = {0,0,0};
     feenox_mesh_subtract(element->node[0]->x, element->node[1]->x, a);
     feenox_mesh_subtract(element->node[0]->x, element->node[2]->x, b);
-    feenox_mesh_normalized_cross(a, b, local_n);
+    feenox_mesh_normalized_cross(a, b, n);
 
   } else if (element->type->dim == 3) {
     feenox_push_error_message("trying to compute the outward normal of a volume (element %d)", element->tag);
@@ -129,33 +129,29 @@ int feenox_mesh_compute_outward_normal(element_t *element, double n[3]) {
   // if there's none (or more than one) then we rely on the element orientation
   if (feenox_mesh_count_element_volumetric_neighbors(element) == 1) {
     // first compute the center of the surface element
-    double surface_center[3];
+    double surface_center[3] = {0,0,0};
     feenox_call(feenox_mesh_compute_element_barycenter(element, surface_center));
 
     // then the center of the volume element
     element_t *volumetric_neighbor = NULL;
     volumetric_neighbor = feenox_mesh_find_element_volumetric_neighbor(element);
 
-    double volumetric_neighbor_center[3];
+    double volumetric_neighbor_center[3] = {0,0,0};
     feenox_call(feenox_mesh_compute_element_barycenter(volumetric_neighbor, volumetric_neighbor_center));
 
     // compute the product between the proposed normal and the difference between these two
     // if the product is positive, invert the normal
-    if (feenox_mesh_subtract_dot(volumetric_neighbor_center, surface_center, local_n) > 0) {
-      local_n[0] = -local_n[0];
-      local_n[1] = -local_n[1];
-      local_n[2] = -local_n[2];
+    if (feenox_mesh_subtract_dot(volumetric_neighbor_center, surface_center, n) > 0) {
+      n[0] = -n[0];
+      n[1] = -n[1];
+      n[2] = -n[2];
     }
   }
 
   // update nx ny and nz
-  feenox_var_value(feenox.mesh.vars.arr_n[0]) = local_n[0];
-  feenox_var_value(feenox.mesh.vars.arr_n[1]) = local_n[1];
-  feenox_var_value(feenox.mesh.vars.arr_n[2]) = local_n[2];
-
-  n[0] = local_n[0];
-  n[1] = local_n[1];
-  n[2] = local_n[2];
+  feenox_var_value(feenox.mesh.vars.arr_n[0]) = n[0];
+  feenox_var_value(feenox.mesh.vars.arr_n[1]) = n[1];
+  feenox_var_value(feenox.mesh.vars.arr_n[2]) = n[2];
 
   return FEENOX_OK;
 }

src/mesh/neighbors.c

@@ -1,7 +1,7 @@
 /*------------ -------------- -------- --- ----- ---   --       -            -
  *  feenox's mesh-related neighbor routines
  *
- *  Copyright (C) 2014--2018 Jeremy Theler
+ *  Copyright (C) 2014--2024 Jeremy Theler
  *
  *  This file is part of feenox.
  *
@@ -61,26 +61,20 @@ element_t *feenox_mesh_find_element_volumetric_neighbor(element_t *this) {
 }
 
 int feenox_mesh_count_element_volumetric_neighbors(element_t *this) {
-
-  // en mallas de primer orden esto sirve para mezclar elementos raros
-  // en segundo hay que hacerlo completo
-//  int target = (this->type->order == 1) ? this->type->dim : this->type->nodes;
-  int target = this->type->nodes;
-
   size_t tags[this->type->faces];
-  for (unsigned int k = 0; k < this->type->faces; k++) {
+  for (int k = 0; k < this->type->faces; k++) {
     tags[k] = 0;
   }
 
   element_ll_t *element_item = NULL;
-  unsigned int index = 0;
+  int index = 0;
   int n = 0;
-  for (unsigned int j = 0; j < this->type->nodes; j++) {
+  for (int j = 0; j < this->type->nodes; j++) {
     LL_FOREACH(this->node[j]->element_list, element_item) {
       if (this->type->dim == (element_item->element->type->dim-1)) {  // los vecinos volumetricos
-        if (mesh_count_common_nodes(this, element_item->element, NULL) >= target) {
+        if (mesh_count_common_nodes(this, element_item->element, NULL) >= this->type->nodes) {
           int exists = 0;
-          for (unsigned int k = 0; exists == 0 && k < index; k++) {
+          for (int k = 0; exists == 0 && k < index; k++) {
             exists |= (element_item->element->tag == tags[k]);
           }
           if (exists == 0) {