Merge branch 'main' into script_update

wortiz · Dec 8, 2023 · 82bf798 · 82bf798
2 parents 82f273b + a02884a
commit 82bf798
Show file tree

Hide file tree

Showing 25 changed files with 215 additions and 116 deletions.
diff --git a/.github/workflows/clang-format-check.yml b/.github/workflows/clang-format-check.yml
@@ -13,7 +13,7 @@ jobs:
     steps:
     - uses: actions/checkout@v3
     - name: Run clang-format style check for C/C++
-      uses: jidicula/clang-format-action@v4.5.0
+      uses: jidicula/clang-format-action@v4.11.0
       with:
-        clang-format-version: '14'
+        clang-format-version: '17'
         check-path: ${{ matrix.path }}
diff --git a/docs/problem_description_file/boundary_conditions/mesh/ca_edge.rst b/docs/problem_description_file/boundary_conditions/mesh/ca_edge.rst
@@ -10,6 +10,11 @@
 **Description / Usage**
 -----------------------
 
+.. warning::
+
+  The angle in this BC is specified in degrees, not radians. This is 
+  different than the 2D equivalent.
+
 **(PCC-EDGE/ROTATED MESH)**
 
 This boundary condition card specifies a constant contact angle on the edge defined by

diff --git a/docs/problem_description_file/boundary_conditions/mesh/ca_edge_int.rst b/docs/problem_description_file/boundary_conditions/mesh/ca_edge_int.rst
@@ -10,6 +10,11 @@
 **Description / Usage**
 -----------------------
 
+.. warning::
+
+  The angle in this BC is specified in degrees, not radians. This is 
+  different than the 2D equivalent.
+
 **(SIC-EDGE/ROTATED MESH)**
 
 This boundary condition card specifies a constant contact angle on the edge defined by

diff --git a/docs/problem_description_file/boundary_conditions/mesh/var_ca_edge.rst b/docs/problem_description_file/boundary_conditions/mesh/var_ca_edge.rst
@@ -10,6 +10,11 @@
 **Description / Usage**
 -----------------------
 
+.. warning::
+
+  The angle in this BC is specified in degrees, not radians. This is 
+  different than the 2D equivalent.
+
 **(SIC-EDGE/ROTATED MESH)**
 
 This card is used to set a variable contact angle on a dynamic three-dimensional contact

diff --git a/include/rf_fem_const.h b/include/rf_fem_const.h
@@ -990,6 +990,7 @@
 #define AC_POSITION      11 // Positional constraint
 #define AC_FLUX_MAT      12
 #define AC_ANGLE         13 // Angle constraint - 2D onl
+#define AC_POSITION_MT   14 // Positional constraint met using material prop.
 
 /* Post Processing Volumetric Integration types - AMC Originally Aug 2013 */
 

diff --git a/src/ac_conti.c b/src/ac_conti.c
@@ -870,7 +870,7 @@ void continue_problem(Comm_Ex *cx, /* array of communications structures */
             evol_local = augc[iAC].evol;
 #ifdef PARALLEL
             if (Num_Proc > 1 && (augc[iAC].Type == AC_VOLUME || augc[iAC].Type == AC_POSITION ||
-                                 augc[iAC].Type == AC_ANGLE)) {
+                                 augc[iAC].Type == AC_ANGLE || augc[iAC].Type == AC_POSITION_MT)) {
               MPI_Allreduce(&evol_local, &evol_global, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
               evol_local = evol_global;
             }
@@ -884,7 +884,7 @@ void continue_problem(Comm_Ex *cx, /* array of communications structures */
             } else if (augc[iAC].Type == AC_VOLUME) {
               DPRINTF(stdout, "\tMT[%4d] VC[%4d]=%10.6e Param=%10.6e\n", augc[iAC].MTID,
                       augc[iAC].VOLID, evol_local, x_AC[iAC]);
-            } else if (augc[iAC].Type == AC_POSITION) {
+            } else if (augc[iAC].Type == AC_POSITION || augc[iAC].Type == AC_POSITION_MT) {
               DPRINTF(stdout, "\tMT[%4d] XY[%4d]=%10.6e Param=%10.6e\n", augc[iAC].MTID,
                       augc[iAC].VOLID, evol_local, x_AC[iAC]);
             } else if (augc[iAC].Type == AC_ANGLE) {

diff --git a/src/ac_hunt.c b/src/ac_hunt.c
@@ -862,7 +862,7 @@ void hunt_problem(Comm_Ex *cx, /* array of communications structures */
             evol_local = augc[iAC].evol;
 #ifdef PARALLEL
             if (Num_Proc > 1 && (augc[iAC].Type == AC_VOLUME || augc[iAC].Type == AC_POSITION ||
-                                 augc[iAC].Type == AC_ANGLE)) {
+                                 augc[iAC].Type == AC_ANGLE || augc[iAC].Type == AC_POSITION_MT)) {
               MPI_Allreduce(&evol_local, &evol_global, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
               evol_local = evol_global;
             }
@@ -876,7 +876,7 @@ void hunt_problem(Comm_Ex *cx, /* array of communications structures */
             } else if (augc[iAC].Type == AC_VOLUME) {
               DPRINTF(stdout, "\tMT[%4d] VC[%4d]=%10.6e Param=%10.6e\n", augc[iAC].MTID,
                       augc[iAC].VOLID, evol_local, x_AC[iAC]);
-            } else if (augc[iAC].Type == AC_POSITION) {
+            } else if (augc[iAC].Type == AC_POSITION || augc[iAC].Type == AC_POSITION_MT) {
               DPRINTF(stdout, "\tMT[%4d] XY[%4d]=%10.6e Param=%10.6e\n", augc[iAC].MTID,
                       augc[iAC].VOLID, evol_local, x_AC[iAC]);
             } else if (augc[iAC].Type == AC_ANGLE) {

diff --git a/src/ac_loca_interface.c b/src/ac_loca_interface.c
@@ -1375,15 +1375,15 @@ int nonlinear_solver_conwrap(double *x, void *con_ptr, int step_num, double lamb
 #endif
           DPRINTF(stdout, "\tMT[%4d] VC[%4d]=%10.6e Param=%10.6e\n", augc[iAC].MTID,
                   augc[iAC].VOLID, evol_local, passdown.x_AC[iAC]);
-        } else if (augc[iAC].Type == AC_POSITION) {
+        } else if (augc[iAC].Type == AC_POSITION || augc[iAC].Type == AC_POSITION_MT) {
           evol_local = augc[iAC].evol;
 #ifdef PARALLEL
           if (Num_Proc > 1) {
             MPI_Allreduce(&evol_local, &evol_global, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
           }
           evol_local = evol_global;
 #endif
-          DPRINTF(stdout, "\tMT[%4d] XY[%4d]=%10.6e Param=%10.6e\n", augc[iAC].MTID,
+          DPRINTF(stdout, "\tMT[%4d] XY[%4d]=%10.6e Param=%13.9e\n", augc[iAC].MTID,
                   augc[iAC].VOLID, evol_local, passdown.x_AC[iAC]);
         } else if (augc[iAC].Type == AC_FLUX) {
           DPRINTF(stdout, "\tBC[%4d] DF[%4d]=%10.6e\n", augc[iAC].BCID, augc[iAC].DFID,

diff --git a/src/ac_particles.c b/src/ac_particles.c
@@ -60,44 +60,6 @@
 /* GOMA include files */
 #define GOMA_AC_PARTICLES_C
 
-/* Global variables extern declared in ac_particles.h. */
-int Particle_Dynamics;                /* global toggle indicating particles are present. */
-enum Particle_Model_t Particle_Model; /* What flavor of particle<->continuum stuff... */
-dbl Particle_Model_Data[MAX_PARTICLE_MODEL_DATA_VALUES]; /* Real values for this model. */
-int Particle_Number;                                     /* number of discrete particles. */
-particle_filename_s Particle_Restart_Filename;           /* restart filename */
-int Particle_Output_Stride;    /* How often to output particle information. */
-dbl Particle_Output_Time_Step; /* Output every these units. */
-int Particle_Max_Time_Steps;   /* max number of particle time steps if steady solution. */
-enum Particle_Output_Format_t Particle_Output_Format; /* What kind of output file? */
-int Particle_Full_Output_Stride;                   /* > 0 => full output every that many steps. */
-particle_filename_s Particle_Full_Output_Filename; /* where to put them. */
-int Particle_Number_Sample_Types;                  /* How many datasets to output? */
-int *Particle_Number_Samples_Existing;             /* How many are tagged for each sample type?*/
-int *Particle_Number_Samples;          /* How many particles to output for dataset #n? */
-int *Particle_Number_Output_Variables; /* How many output vars for each sample. */
-particle_variable_s *
-    *Particle_Output_Variables; /* List of variable indices to output for dataset #n */
-particle_filename_s *Particle_Filename_Template; /* Template of where to put the data... */
-
-dbl Particle_Density;         /* Density of particle in problem units */
-dbl Particle_Radius;          /* Radius of particle in problem units. */
-dbl Particle_Ratio;           /* Real/computational particle ratio. */
-int Particle_Show_Debug_Info; /* Show particle debug info. */
-enum Particle_Domain_t Particle_Creation_Domain;
-enum Particle_Domain_t Particle_Move_Domain;
-particle_filename_s Particle_Creation_Domain_Filename;
-particle_s Particle_Creation_Domain_Name;
-particle_filename_s Particle_Move_Domain_Filename;
-particle_s Particle_Move_Domain_Name;
-dbl Particle_Creation_Domain_Reals[MAX_DOMAIN_REAL_VALUES];
-dbl Particle_Move_Domain_Reals[MAX_DOMAIN_REAL_VALUES];
-dbl xi_boundary_tolerances[3] = {XI_BOUNDARY_TOLERANCE0, XI_BOUNDARY_TOLERANCE1,
-                                 XI_BOUNDARY_TOLERANCE2};
-
-int Particle_Number_PBCs; /* Number of particle-related sideset BC's. */
-PBC_t *PBCs;              /* Particle boundary condition structures. */
-
 /* Global variables that reside entirely within this file. */
 static particle_t *particles_to_do, *particles_to_send;
 static particle_t **element_particle_list_head;
@@ -128,8 +90,6 @@ static dbl *static_xdot;
 static dbl *static_xdot_old;
 static dbl *static_resid_vector;
 
-element_particle_info_t *element_particle_info;
-
 static dbl my_volume;
 static dbl *el_volume;
 static FILE **pa_fp;

diff --git a/src/bc/rotate_util.c b/src/bc/rotate_util.c
@@ -294,7 +294,7 @@ goma_corner_coordinate_system(goma_normal **normals, int n_normals, goma_normal
 
   return GOMA_ERROR;
 
-corner_coord_critical_found : {
+corner_coord_critical_found: {
   // find closest to x for first
   int first_crit = -1;
   double max_dot = 0;

diff --git a/src/dp_ghost.cpp b/src/dp_ghost.cpp
@@ -16,7 +16,6 @@
 #include <cstdio>
 #include <cstring>
 #include <exodusII.h>
-#include <ext/alloc_traits.h>
 #include <iterator>
 #include <memory>
 #include <mpi.h>

diff --git a/src/globals.c b/src/globals.c
@@ -13,6 +13,7 @@
 \************************************************************************/
 /* File containing the many of Goma's global variables, moved for -fno-common */
 
+#include "ac_particles.h"
 #include "rf_fem.h"
 #include "rf_fem_const.h"
 #include "rf_node_const.h"
@@ -297,3 +298,46 @@ struct Boundary_Condition *BC_Types;
 struct Rotation_Specs *ROT_Types;
 
 NODE_INFO_STRUCT **Nodes = NULL;
+
+struct elem_side_bc_struct ***First_Elem_Side_BC_Array;
+struct elem_edge_bc_struct ***First_Elem_Edge_BC_Array;
+
+/* Global variables extern declared in ac_particles.h. */
+int Particle_Dynamics;                /* global toggle indicating particles are present. */
+enum Particle_Model_t Particle_Model; /* What flavor of particle<->continuum stuff... */
+dbl Particle_Model_Data[MAX_PARTICLE_MODEL_DATA_VALUES]; /* Real values for this model. */
+int Particle_Number;                                     /* number of discrete particles. */
+particle_filename_s Particle_Restart_Filename;           /* restart filename */
+int Particle_Output_Stride;    /* How often to output particle information. */
+dbl Particle_Output_Time_Step; /* Output every these units. */
+int Particle_Max_Time_Steps;   /* max number of particle time steps if steady solution. */
+enum Particle_Output_Format_t Particle_Output_Format; /* What kind of output file? */
+int Particle_Full_Output_Stride;                   /* > 0 => full output every that many steps. */
+particle_filename_s Particle_Full_Output_Filename; /* where to put them. */
+int Particle_Number_Sample_Types;                  /* How many datasets to output? */
+int *Particle_Number_Samples_Existing;             /* How many are tagged for each sample type?*/
+int *Particle_Number_Samples;          /* How many particles to output for dataset #n? */
+int *Particle_Number_Output_Variables; /* How many output vars for each sample. */
+particle_variable_s *
+    *Particle_Output_Variables; /* List of variable indices to output for dataset #n */
+particle_filename_s *Particle_Filename_Template; /* Template of where to put the data... */
+
+dbl Particle_Density;         /* Density of particle in problem units */
+dbl Particle_Radius;          /* Radius of particle in problem units. */
+dbl Particle_Ratio;           /* Real/computational particle ratio. */
+int Particle_Show_Debug_Info; /* Show particle debug info. */
+enum Particle_Domain_t Particle_Creation_Domain;
+enum Particle_Domain_t Particle_Move_Domain;
+particle_filename_s Particle_Creation_Domain_Filename;
+particle_s Particle_Creation_Domain_Name;
+particle_filename_s Particle_Move_Domain_Filename;
+particle_s Particle_Move_Domain_Name;
+dbl Particle_Creation_Domain_Reals[MAX_DOMAIN_REAL_VALUES];
+dbl Particle_Move_Domain_Reals[MAX_DOMAIN_REAL_VALUES];
+dbl xi_boundary_tolerances[3] = {XI_BOUNDARY_TOLERANCE0, XI_BOUNDARY_TOLERANCE1,
+                                 XI_BOUNDARY_TOLERANCE2};
+
+int Particle_Number_PBCs; /* Number of particle-related sideset BC's. */
+PBC_t *PBCs;              /* Particle boundary condition structures. */
+
+element_particle_info_t *element_particle_info;
diff --git a/src/md_timer.c b/src/md_timer.c
@@ -67,6 +67,11 @@
 #include <sys/types.h>
 #endif
 
+#ifdef __MACH__
+#include <sys/times.h>
+#include <time.h>
+#endif
+
 #define GOMA_MD_TIMER_C
 #include "md_timer.h"
 #include "std.h"

diff --git a/src/mm_augc_util.c b/src/mm_augc_util.c
@@ -128,7 +128,8 @@ void load_extra_unknownsAC(int iAC,     /* ID NUMBER OF AC'S */
   }
 
   if (augc[iAC].Type == AC_USERBC || augc[iAC].Type == AC_VOLUME || augc[iAC].Type == AC_FLUX ||
-      augc[iAC].Type == AC_LS_VEL || augc[iAC].Type == AC_POSITION || augc[iAC].Type == AC_ANGLE) {
+      augc[iAC].Type == AC_LS_VEL || augc[iAC].Type == AC_POSITION || augc[iAC].Type == AC_ANGLE ||
+      augc[iAC].Type == AC_POSITION_MT) {
 
     ibc = augc[iAC].BCID;
     idf = augc[iAC].DFID;
@@ -1962,7 +1963,7 @@ int std_aug_cond(int iAC,
     // Formulate and store the residual for the augmented condition
     gAC[iAC] = (inventory - augc[iAC].CONSTV);
 
-  } else if (augc[iAC].Type == AC_POSITION) {
+  } else if (augc[iAC].Type == AC_POSITION || augc[iAC].Type == AC_POSITION_MT) {
     inventory = getPositionAC(augc + iAC, cAC[iAC], mf_args->x, mf_args->exo);
 #ifdef PARALLEL
     if (Num_Proc > 1) {
@@ -2571,6 +2572,23 @@ static int estimate_dAC_ALC(
       }
       break;
 
+    case AC_POSITION:
+      // here you are supposed to calculate the residual at the + perturbation
+      // I believe this can be done for the position type ac by copying the structure above (line
+      // 1992) this should calculate the residual at the updated x
+
+      inventory = getPositionAC(augc + jAC, dAC[jAC], mf_args->x, mf_args->exo);
+#ifdef PARALLEL
+      if (Num_Proc > 1) {
+        MPI_Allreduce(&inventory, &global_inventory, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        inventory = global_inventory;
+      }
+#endif
+      // Formulate and store the perturbed residual for the augmented condition
+      // seems correct. -BW 04/18/2022 + WO 04/21/2023
+      res_p[jAC] = inventory - augc[jAC].CONSTV;
+      break;
+
     case AC_USERBC:
     case AC_USERMAT:
       load_extra_unknownsAC(jAC, x_AC, cx, mf_args->exo, mf_args->dpi);
@@ -2661,6 +2679,17 @@ static int estimate_dAC_ALC(
         res_m[jAC] = inventory - (-BC_Types[ibc].BC_Data_Float[idf]);
       }
       break;
+    case AC_POSITION:
+      // just do it again
+      inventory = getPositionAC(augc + jAC, dAC[jAC], mf_args->x, mf_args->exo);
+#ifdef PARALLEL
+      if (Num_Proc > 1) {
+        MPI_Allreduce(&inventory, &global_inventory, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        inventory = global_inventory;
+      }
+#endif
+      res_m[jAC] = inventory - augc[jAC].CONSTV;
+      break;
 
     case AC_USERBC:
     case AC_USERMAT:

diff --git a/src/mm_fill.c b/src/mm_fill.c
@@ -110,6 +110,9 @@
 
 #define GOMA_MM_FILL_C
 
+extern struct elem_side_bc_struct ***First_Elem_Side_BC_Array;
+extern struct elem_edge_bc_struct ***First_Elem_Edge_BC_Array;
+
 /*
  * Global variables defined here. Declared frequently via rf_bc.h
  */
@@ -5453,7 +5456,7 @@ int checkfinite(const char *file, const int line, const char *message) {
     peqn = upd->ep[pg->imtrx][eqn];
     if (peqn != -1) {
       for (i = 0; i < ei[pg->imtrx]->dof[eqn]; i++) {
-        j = finite(lec->R[LEC_R_INDEX(peqn, i)]);
+        j = isfinite(lec->R[LEC_R_INDEX(peqn, i)]);
         if (!j) {
           fprintf(stderr, "lec->R[%s][edof=%d] = %g\n", EQ_Name[eqn].name1, i,
                   lec->R[LEC_R_INDEX(peqn, i)]);

diff --git a/src/mm_fill_rs.c b/src/mm_fill_rs.c
@@ -1967,41 +1967,41 @@ void f_kinematic_displacement_bc(double func[DIM],
     }
 #endif
 
-  if (af->Assemble_Jacobian) {
-    for (kdir = 0; kdir < pd->Num_Dim; kdir++) {
-      for (p = 0; p < pd->Num_Dim; p++) {
-        var = MESH_DISPLACEMENT1 + p;
-        if (pd->v[pg->imtrx][var]) {
-          for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) {
-            phi_j = bf[var]->phi[j];
-            d_func[0][var][j] += ((fv->d[kdir] - base_displacement[kdir]) -
-                                  (fv->d_rs[kdir] - base_displacement_rs[kdir])) *
-                                     fv->dsnormal_dx[kdir][p][j] +
-                                 (delta(kdir, p) - dns_dX[kdir][p]) * phi_j * fv->snormal[kdir];
+    if (af->Assemble_Jacobian) {
+      for (kdir = 0; kdir < pd->Num_Dim; kdir++) {
+        for (p = 0; p < pd->Num_Dim; p++) {
+          var = MESH_DISPLACEMENT1 + p;
+          if (pd->v[pg->imtrx][var]) {
+            for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) {
+              phi_j = bf[var]->phi[j];
+              d_func[0][var][j] += ((fv->d[kdir] - base_displacement[kdir]) -
+                                    (fv->d_rs[kdir] - base_displacement_rs[kdir])) *
+                                       fv->dsnormal_dx[kdir][p][j] +
+                                   (delta(kdir, p) - dns_dX[kdir][p]) * phi_j * fv->snormal[kdir];
+            }
           }
         }
-      }
 
-      for (p = 0; p < pd->Num_Dim; p++) {
-        var = SOLID_DISPLACEMENT1 + p;
-        if (pd->v[pg->imtrx][var]) {
-          for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) {
-            phi_j = bf[var]->phi[j];
-            d_func[0][var][j] -= (delta(kdir, p) + dns_drs[kdir][p]) * phi_j * fv->snormal[kdir];
+        for (p = 0; p < pd->Num_Dim; p++) {
+          var = SOLID_DISPLACEMENT1 + p;
+          if (pd->v[pg->imtrx][var]) {
+            for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) {
+              phi_j = bf[var]->phi[j];
+              d_func[0][var][j] -= (delta(kdir, p) + dns_drs[kdir][p]) * phi_j * fv->snormal[kdir];
+            }
           }
         }
       }
-    }
-  } /* end of if Assemble_Jacobian */
+    } /* end of if Assemble_Jacobian */
 
-  /* Calculate the residual contribution	*/
+    /* Calculate the residual contribution	*/
 
-  for (kdir = 0; kdir < pd->Num_Dim; kdir++) {
-    *func +=
-        ((fv->d[kdir] - base_displacement[kdir]) - (fv->d_rs[kdir] - base_displacement_rs[kdir])) *
-        fv->snormal[kdir];
+    for (kdir = 0; kdir < pd->Num_Dim; kdir++) {
+      *func += ((fv->d[kdir] - base_displacement[kdir]) -
+                (fv->d_rs[kdir] - base_displacement_rs[kdir])) *
+               fv->snormal[kdir];
+    }
   }
-}
 
 } /* END of routine f_kinematic_displacement_bc  */