Add functionality for volumetric coupling of material models

Makefile

@@ -52,7 +52,7 @@ LIBS = \
 #CFLAGS = -g -Wall -std=c++11 -O0 -fopenmp $(INCLUDES) -DARCH="Linux" -DSPOOLES -DARPACK -DMATRIXSTORAGE
 #FFLAGS = -g -Wall -O0 -fopenmp $(INCLUDES)
 
-CFLAGS = -Wall -O3 -fopenmp $(INCLUDES) -DARCH="Linux" -DSPOOLES -DARPACK -DMATRIXSTORAGE -DUSE_MT
+CFLAGS = -Wall -O0 -fopenmp $(INCLUDES) -DARCH="Linux" -DSPOOLES -DARPACK -DMATRIXSTORAGE -DUSE_MT -fbounds-check -g
 
 # OS-specific options
 UNAME_S := $(shell uname -s)
@@ -62,7 +62,7 @@ else
 	CC = mpicc
 endif
 
-FFLAGS = -Wall -O3 -fopenmp $(INCLUDES) ${ADDITIONAL_FFLAGS}
+FFLAGS = -Wall -O0 -fopenmp $(INCLUDES) ${ADDITIONAL_FFLAGS} -fallow-argument-mismatch -fbounds-check -g
 # Note for GCC 10 or newer: add -fallow-argument-mismatch in the above flags
 FC = mpifort
 # FC = mpif90
@@ -73,10 +73,8 @@ include $(CCX)/Makefile.inc
 SCCXMAIN = ccx_$(CCX_VERSION).c
 
 # Append additional sources
-SCCXC += nonlingeo_precice.c dyna_precice.c CCXHelpers.c PreciceInterface.c
-SCCXF += getflux.f getkdeltatemp.f
-
-
+SCCXC += nonlingeo_precice.c dyna_precice.c CCXHelpers.c PreciceInterface.c linstatic_precice.c
+SCCXF += getflux.f getkdeltatemp.f getc3d8elementgausspointcoords.f multiscale_routines.f
 
 # Source files in this folder and in the adapter directory
 $(OBJDIR)/%.o : %.c
@@ -104,7 +102,7 @@ OCCXC += $(OBJDIR)/ConfigReader.o $(OBJDIR)/2D3DCoupling.o $(OBJDIR)/OutputBuffe
 
 
 $(OBJDIR)/ccx_preCICE: $(OBJDIR) $(OCCXMAIN) $(OBJDIR)/ccx_$(CCX_VERSION).a
-	$(FC) -fopenmp -Wall -O3 -o $@ $(OCCXMAIN) $(OBJDIR)/ccx_$(CCX_VERSION).a $(LIBS)
+	$(FC) -fopenmp -Wall -O0 -g -o $@ $(OCCXMAIN) $(OBJDIR)/ccx_$(CCX_VERSION).a $(LIBS)
 
 $(OBJDIR)/ccx_$(CCX_VERSION).a: $(OCCXF) $(OCCXC)
 	ar vr $@ $?

adapter/CCXHelpers.c

@@ -68,6 +68,17 @@ void getSurfaceElementsAndFaces(ITG setID, ITG *ialset, ITG *istartset, ITG *ien
   }
 }
 
+void getElementsIDs(ITG setID, ITG *ialset, ITG *istartset, ITG *iendset, ITG *elements)
+{
+
+  ITG i, k = 0;
+
+  for (i = istartset[setID] - 1; i < iendset[setID]; i++) {
+    elements[k] = ialset[i];
+    k++;
+  }
+}
+
 void getNodeCoordinates(ITG *nodes, ITG numNodes, int dim, double *co, double *v, int mt, double *coordinates)
 {
 
@@ -265,7 +276,7 @@ void getXloadIndices(char const *loadType, ITG *elementIDs, ITG *faceIDs, ITG nu
   char faceLabel[] = {'x', 'x', '\0'};
 
   /* Face number is prefixed with 'S' if it is DFLUX boundary condition
-   * and with 'F' if it is a FILM boundary condition */
+	 * and with 'F' if it is a FILM boundary condition */
   if (strcmp(loadType, "DFLUX") == 0) {
     faceLabel[0] = (char) 'S';
   } else if (strcmp(loadType, "FILM") == 0) {
@@ -355,22 +366,22 @@ void getXforcIndices(ITG *nodes, ITG numNodes, int nforc, int *ikforc, int *ilfo
   ITG i;
 
   for (i = 0; i < numNodes; i++) {
-    // x-direction
+    //x-direction
     int idof = 8 * (nodes[i] - 1) + 1; // 1 for x force DOF
     int k;
     FORTRAN(nident, (ikforc, &idof, &nforc, &k));
     k -= 1;                              // Adjust because of FORTRAN indices
     int m               = ilforc[k] - 1; // Adjust because of FORTRAN indices
     xforcIndices[3 * i] = m;
 
-    // y-direction
+    //y-direction
     idof = 8 * (nodes[i] - 1) + 2; // 2 for y force DOF
     FORTRAN(nident, (ikforc, &idof, &nforc, &k));
     k -= 1;                                  // Adjust because of FORTRAN indices
     m                       = ilforc[k] - 1; // Adjust because of FORTRAN indices
     xforcIndices[3 * i + 1] = m;
 
-    // z-direction
+    //z-direction
     idof = 8 * (nodes[i] - 1) + 3; // 3 for z force DOF
     FORTRAN(nident, (ikforc, &idof, &nforc, &k));
     k -= 1;                                  // Adjust because of FORTRAN indices
@@ -389,11 +400,11 @@ void getXforcIndices(ITG *nodes, ITG numNodes, int nforc, int *ikforc, int *ilfo
 int getXloadIndexOffset(enum xloadVariable xloadVar)
 {
   /*
-   * xload is the CalculiX array where the DFLUX and FILM boundary conditions are stored
-   * the array has two components:
-   * - the first component corresponds to the flux value and the heat transfer coefficient
-   * - the second component corresponds to the sink temperature
-   * */
+	 * xload is the CalculiX array where the DFLUX and FILM boundary conditions are stored
+	 * the array has two components:
+	 * - the first component corresponds to the flux value and the heat transfer coefficient
+	 * - the second component corresponds to the sink temperature
+	 * */
   switch (xloadVar) {
   case DFLUX:
     return 0;
@@ -409,6 +420,20 @@ int getXloadIndexOffset(enum xloadVariable xloadVar)
   }
 }
 
+void getElementStrain(int strainIdx, int *mi, int nelem, double *eei, double *strainData)
+{
+  int i, count, idx;
+  // Loop through all element and respective gauss points
+  count = 0;
+  for (i = 0; i < mi[0] * nelem; i++) {
+    idx                   = i * 6 + strainIdx;
+    strainData[count]     = eei[idx];
+    strainData[count + 1] = eei[idx + 1];
+    strainData[count + 2] = eei[idx + 2];
+    count                 = count + 3;
+  }
+}
+
 void setXload(double *xload, int *xloadIndices, double *values, int numValues, enum xloadVariable xloadVar)
 {
   ITG i;
@@ -512,6 +537,17 @@ bool isQuasi2D3D(const int quasi2D3D)
   return quasi2D3D == 1;
 }
 
+void setDoubleArrayZero(double *values, const int length, const int dim)
+{
+  ITG i, j;
+
+  for (i = 0; i < length; i++) {
+    for (j = 0; j < dim; j++) {
+      values[i * dim + j] = 0.0;
+    }
+  }
+}
+
 /* Errors messages */
 
 void nodeSetNotFoundError(char const *setName)

adapter/CCXHelpers.h

@@ -51,17 +51,23 @@ enum CouplingDataType { TEMPERATURE,
                         DISPLACEMENTS,
                         DISPLACEMENTDELTAS,
                         VELOCITIES,
-                        POSITIONS };
-
-/**
- * @brief Type of element used for faces mesh, where we assume only one type of element is used.
- *  TETRAHEDRA      - C3D4 or C3D10 element
- *  HEXAHEDRA        - C3D8 or C3D20 elment
- *  INVALID_ELEMENT - Anything else
- */
-enum ElemType { TETRAHEDRA,
-                HEXAHEDRA,
-                INVALID_ELEMENT,
+                        POSITIONS,
+                        MACRO_IP_ID,
+                        RVE_ID,
+                        MOD_ID,
+                        RUC_SIZE,
+                        CONV_FLAG,
+                        STRAIN1TO3,
+                        STRAIN4TO6,
+                        STRESS1TO3,
+                        STRESS4TO6,
+                        CMAT1,
+                        CMAT2,
+                        CMAT3,
+                        CMAT4,
+                        CMAT5,
+                        CMAT6,
+                        CMAT7,
 };
 
 /**
@@ -102,6 +108,16 @@ ITG getNumSetElements(ITG setID, ITG *istartset, ITG *iendset);
  */
 void getSurfaceElementsAndFaces(ITG setID, ITG *ialset, ITG *istartset, ITG *iendset, ITG *elements, ITG *faces);
 
+/**
+ * @brief Gets the element IDs given a set ID
+ * @param setID: input set id
+ * @param ialset: CalculiX variable
+ * @param istartset: CalculiX variable
+ * @param iendset: CalculiX variable
+ * @param elements: output element IDs
+ */
+void getElementsIDs(ITG setID, ITG *ialset, ITG *istartset, ITG *iendset, ITG *elements);
+
 /**
  * @brief Gets the coordinates of a list of input node IDs
  * @param nodes: input node IDs
@@ -241,6 +257,16 @@ void getXbounIndices(ITG *nodes, ITG numNodes, int nboun, int *ikboun, int *ilbo
  */
 void getXforcIndices(ITG *nodes, ITG numNodes, int nforc, int *ikforc, int *ilforc, int *xforcIndices);
 
+/**
+ * @brief Gets the strain values at each Gauss point of each element
+ * @param strainIdx: CalculiX variable for the index of the strain values
+ * @param mi: CalculiX variable for the number of integration points
+ * @param nelem: CalculiX variable for the number of elements
+ * @param eei: CalculiX array for the element integration information
+ * @param strainData: CalculiX array for the strain values
+ */
+void getElementStrain(int strainIdx, int *mi, int nelem, double *eei, double *strainData);
+
 /**
  * @brief Modifies the values of a DFLUX or FILM boundary condition
  * @param xload: CalculiX array for the loads

adapter/ConfigReader.cpp

@@ -53,6 +53,10 @@ void ConfigReader_Read(char const *configFilename, char const *participantName,
       interface.facesMeshName = strdup(config["participants"][participantName]["interfaces"][i]["mesh"].as<std::string>().c_str());
     }
 
+    if (config["participants"][participantName]["interfaces"][i]["elements"]) {
+      interface.elementsMeshName = strdup(config["participants"][participantName]["interfaces"][i]["elements"].as<std::string>().c_str());
+    }
+
     std::string patchName = config["participants"][participantName]["interfaces"][i]["patch"].as<std::string>();
     std::transform(patchName.begin(), patchName.end(), patchName.begin(), toupper);
     interface.patchName = strdup(patchName.c_str());

adapter/ConfigReader.h

@@ -13,6 +13,7 @@
 typedef struct InterfaceConfig {
   char * facesMeshName;
   char * nodesMeshName;
+  char * elementsMeshName;
   char * patchName;
   int    map;
   int    numWriteData;