Functionality for volumetric coupling on quadrature points of a mesh

Makefile

@@ -77,7 +77,7 @@ SCCXMAIN = ccx_$(CCX_VERSION).c
 
 # Append additional sources
 SCCXC += nonlingeo_precice.c dyna_precice.c CCXHelpers.c PreciceInterface.c
-SCCXF += getflux.f getkdeltatemp.f
+SCCXF += getflux.f getkdeltatemp.f getc3d8elementgausspointcoords.f getc3d4elementgausspointcoords.f
 
 
 

adapter/CCXHelpers.c

@@ -68,6 +68,16 @@ 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)
 {
 
@@ -219,15 +229,6 @@ void getHexaFaceCenters(ITG *elements, ITG *faces, ITG numElements, ITG *kon, IT
   }
 }
 
-/*
-   void getSurfaceGaussPoints(int setID, ITG * co, ITG istartset, ITG iendset, ITG * ipkon, ITG * lakon, ITG * kon, ITG * ialset, double * coords) {
-
-        int iset = setID + 1; // plus one because of fortran indices
-        FORTRAN(getgausspointscoords,(co,&iset,istartset,iendset,ipkon,lakon,kon,ialset, coords));
-
-   }
- */
-
 void getTetraFaceNodes(ITG *elements, ITG *faces, ITG *nodes, ITG numElements, ITG numNodes, ITG *kon, ITG *ipkon, int *tetraFaceNodes)
 {
   // Assume all tetra elements -- maybe implement checking later...

adapter/CCXHelpers.h

@@ -91,6 +91,16 @@ ITG getSetID(char const *setName, char const *set, ITG nset);
  */
 ITG getNumSetElements(ITG setID, ITG *istartset, ITG *iendset);
 
+/**
+ * @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 element and face IDs given a set ID
  * @param setID: input set id

adapter/ConfigReader.cpp

@@ -55,6 +55,10 @@ void ConfigReader_Read(char const *configFilename, char const *participantName,
       interface.facesMeshName = NULL;
     }
 
+    if (config["participants"][participantName]["interfaces"][i].contains("elements-mesh")) {
+      interface.elementsMeshName = strdup(config["participants"][participantName]["interfaces"][i]["elements-mesh"].get_value<std::string>().c_str());
+    }
+
     if (config["participants"][participantName]["interfaces"][i].contains("mesh")) {
       interface.facesMeshName = strdup(config["participants"][participantName]["interfaces"][i]["mesh"].get_value<std::string>().c_str());
     }

adapter/ConfigReader.h

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

adapter/PreciceInterface.c

@@ -450,8 +450,8 @@ void Precice_FreeData(SimulationData *sim)
 void PreciceInterface_Create(PreciceInterface *interface, SimulationData *sim, InterfaceConfig const *config)
 {
   // Deduce configured dimensions
-  if (config->nodesMeshName == NULL && config->facesMeshName == NULL) {
-    printf("ERROR: You need to define either a face or a nodes mesh. Check the adapter configuration file.\n");
+  if (config->nodesMeshName == NULL && config->facesMeshName == NULL && config->elementsMeshName == NULL) {
+    printf("ERROR: You need to define either a face, nodes, or elements mesh. Check the adapter configuration file.\n");
     exit(EXIT_FAILURE);
   }
   if (config->nodesMeshName && config->facesMeshName) {
@@ -473,6 +473,8 @@ void PreciceInterface_Create(PreciceInterface *interface, SimulationData *sim, I
 
   // Initialize pointers as NULL
   interface->elementIDs            = NULL;
+  interface->elemIPID              = NULL;
+  interface->elemIPCoordinates     = NULL;
   interface->faceIDs               = NULL;
   interface->faceCenterCoordinates = NULL;
   interface->preciceFaceCenterIDs  = NULL;
@@ -545,6 +547,14 @@ void PreciceInterface_Create(PreciceInterface *interface, SimulationData *sim, I
     interface->couplingMeshName = interface->faceCentersMeshName;
   }
 
+  // Element mesh
+  interface->elementsMeshName = NULL;
+  if (config->elementsMeshName) {
+    interface->elementsMeshName = strdup(config->elementsMeshName);
+    PreciceInterface_ConfigureElementsMesh(interface, sim);
+    interface->couplingMeshName = interface->elementsMeshName;
+  }
+
   PreciceInterface_ConfigureCouplingData(interface, sim, config);
 }
 
@@ -567,6 +577,48 @@ static enum ElemType findSimulationMeshType(SimulationData *sim)
   return INVALID_ELEMENT;
 }
 
+void PreciceInterface_ConfigureElementsMesh(PreciceInterface *interface, SimulationData *sim)
+{
+  char *elementSetName    = interface->name;
+  interface->elementSetID = getSetID(elementSetName, sim->set, sim->nset);
+  interface->numElements  = getNumSetElements(interface->elementSetID, sim->istartset, sim->iendset);
+
+  interface->elementIDs = malloc(interface->numElements * sizeof(ITG));
+  getElementsIDs(interface->elementSetID, sim->ialset, sim->istartset, sim->iendset, interface->elementIDs);
+
+  interface->numIPTotal        = sim->mi[0] * interface->numElements; // Number of Gauss points per element * number of elements
+  interface->elemIPCoordinates = malloc(interface->numIPTotal * 3 * sizeof(double));
+
+  interface->elemIPID = malloc(interface->numIPTotal * sizeof(int));
+  for (int j = 0; j < interface->numIPTotal; j++) {
+    interface->elemIPID[j] = j;
+  }
+
+  int numElements = interface->numElements;
+
+  enum ElemType elemType = findSimulationMeshType(sim);
+
+  if (elemType == TETRAHEDRA) {
+    FORTRAN(getc3d4elementgausspointcoords, (&numElements,
+                                             interface->elementIDs,
+                                             sim->co,
+                                             sim->kon,
+                                             sim->ipkon,
+                                             interface->elemIPCoordinates));
+  } else if (elemType == HEXAHEDRA) {
+    FORTRAN(getc3d8elementgausspointcoords, (&numElements,
+                                             interface->elementIDs,
+                                             sim->co,
+                                             sim->kon,
+                                             sim->ipkon,
+                                             interface->elemIPCoordinates));
+  } else {
+    supportedElementError();
+  }
+
+  precicec_setMeshVertices(interface->elementsMeshName, interface->numIPTotal, interface->elemIPCoordinates, interface->elemIPID);
+}
+
 void PreciceInterface_ConfigureFaceCentersMesh(PreciceInterface *interface, SimulationData *sim)
 {
   // printf("Entering ConfigureFaceCentersMesh \n");
@@ -803,7 +855,10 @@ void PreciceInterface_FreeData(PreciceInterface *preciceInterface)
   free(preciceInterface->readData);
   free(preciceInterface->writeData);
   free(preciceInterface->elementIDs);
+  free(preciceInterface->elemIPID);
+  free(preciceInterface->elemIPCoordinates);
   free(preciceInterface->faceIDs);
+  free(preciceInterface->nodeIDs);
   free(preciceInterface->preciceFaceCenterIDs);
   free(preciceInterface->faceCenterCoordinates);
   free(preciceInterface->nodeCoordinates);
@@ -822,6 +877,7 @@ void PreciceInterface_FreeData(PreciceInterface *preciceInterface)
   // Mesh names
   free(preciceInterface->faceCentersMeshName);
   free(preciceInterface->nodesMeshName);
+  free(preciceInterface->elementsMeshName);
 
   // Data names
   free(preciceInterface->displacementDeltas);

adapter/PreciceInterface.h

@@ -33,10 +33,15 @@ typedef struct PreciceInterface {
   int          nodeSetID;
   int *        preciceNodeIDs;
   char *       nodesMeshName;
+  char *       elementsMeshName;
 
-  // Interface face elements
   int     numElements;
   int *   elementIDs;
+  int *   elemIPID;
+  int     elementSetID;
+  double *elemIPCoordinates;
+  int     numIPTotal;
+
   int *   faceIDs;
   double *faceCenterCoordinates;
   int     faceSetID;
@@ -277,6 +282,13 @@ void sendFaceCentersVertices(PreciceInterface *interface);
  */
 void PreciceInterface_ConfigureNodesMesh(PreciceInterface *interface, SimulationData *sim);
 
+/**
+ * @brief Configures the elements mesh and calls setMeshVertices on preCICE
+ * @param interface
+ * @param sim: Structure with CalculiX data
+ */
+void PreciceInterface_ConfigureElementsMesh(PreciceInterface *interface, SimulationData *sim);
+
 /**
  * @brief Terminate execution if the nodes mesh ID is not valid
  * @param interface

getc3d4elementgausspointcoords.f

@@ -0,0 +1,86 @@
+!
+!     CalculiX - A 3-dimensional finite element program
+!              Copyright (C) 1998-2015 Guido Dhondt
+!
+!     This program is free software; you can redistribute it and/or
+!     modify it under the terms of the GNU General Public License as
+!     published by the Free Software Foundation(version 2);
+!     
+!
+!     This program is distributed in the hope that it will be useful,
+!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
+!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
+!     GNU General Public License for more details.
+!
+!     You should have received a copy of the GNU General Public License
+!     along with this program; if not, write to the Free Software
+!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+!
+!     This function is a copy of the function for calculation and printout of the lift and drag forces
+!
+      subroutine getc3d4elementgausspointcoords(nelem, elem_ids,
+     &    co, kon, ipkon, elem_gp_coord)
+
+         implicit none
+
+         !Input/Output variables
+         integer     :: nelem            ! Number of elements
+         integer     :: elem_ids(*)      ! Element IDs, need to be adjusted for Fortran
+         real(8)     :: co(3,*)          ! Nodal coordinates of all nodes
+         integer     :: kon(*)
+         integer     :: ipkon(*)
+
+         real(8)     :: elem_gp_coord(*) ! Gauss point coords
+
+         ! Internal variables
+         integer(4)     :: iel, el_id, indexe, i, j, k, l
+         integer(4)     :: iflag, idx, konl(4), gp_id
+         REAL(8)        :: xl(3,4), xi, et, ze, xsj, shp(4,20)
+
+         include "gauss.f"
+
+         data iflag /1/
+
+         ! Initialize gauss point coordinates to zero
+         do l = 1, nelem*3
+            elem_gp_coord(l) = 0.0
+         end do
+
+         do iel = 1, nelem
+            el_id = elem_ids(iel)
+
+            indexe=ipkon(el_id)
+
+            ! Local nodal coordinates
+            do i = 1, 4
+               konl(i) = kon(indexe+i)
+               do j = 1, 3
+                  xl(j,i)=co(j,konl(i))
+               end do
+            end do
+
+            ! One Gauss point per element
+            gp_id = 1
+
+            ! gauss3d4: tet, 1 integration point
+            xi = gauss3d4(1,gp_id)
+            et = gauss3d4(2,gp_id)
+            ze = gauss3d4(3,gp_id)
+
+            ! Get the shape function
+            call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
+
+            ! Calculate the Gauss point coordinates
+            do k = 1, 3
+               idx = (iel - 1) * 3 + k
+               do l = 1, 4
+                  elem_gp_coord(idx)=elem_gp_coord(idx)
+     &                  +xl(k,l)*shp(4,l)
+               end do
+            end do
+
+         end do
+
+         return
+
+      end subroutine getc3d4elementgausspointcoords

getc3d8elementgausspointcoords.f

@@ -0,0 +1,91 @@
+!
+!     CalculiX - A 3-dimensional finite element program
+!              Copyright (C) 1998-2015 Guido Dhondt
+!
+!     This program is free software; you can redistribute it and/or
+!     modify it under the terms of the GNU General Public License as
+!     published by the Free Software Foundation(version 2);
+!     
+!
+!     This program is distributed in the hope that it will be useful,
+!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
+!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
+!     GNU General Public License for more details.
+!
+!     You should have received a copy of the GNU General Public License
+!     along with this program; if not, write to the Free Software
+!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+!
+!     This function is a copy of the function for calculation and printout of the lift and drag forces
+!
+      subroutine getc3d8elementgausspointcoords(nelem, elem_ids,
+     &    co, kon, ipkon, elem_gp_coord)
+
+         implicit none
+
+         !Input/Output variables
+         integer     :: nelem            ! Number of elements
+         integer     :: elem_ids(*)      ! Element IDs, need to be adjusted for Fortran
+         real(8)     :: co(3,*)          ! Nodal coordinates of all nodes
+         integer     :: kon(*)
+         integer     :: ipkon(*)
+
+         real(8)     :: elem_gp_coord(*) ! GP coords
+
+         ! Internal variables
+         integer(4)     :: iel, el_id, indexe, i, j, k, l
+         integer(4)     :: iflag, idx, konl(8), gp_id
+         real(8)        :: xl(3,8), xi, et, ze, xsj, shp(4,20)
+
+         include "gauss.f"
+
+         data iflag /1/
+
+         ! Initialize gauss point coordinates to zero
+         do l = 1, nelem*8*3
+            elem_gp_coord(l) = 0.0
+         end do
+
+         do iel = 1, nelem
+            el_id = elem_ids(iel)
+
+            indexe=ipkon(el_id)
+
+            ! connectivity
+            do i = 1, 8
+               konl(i)=kon(indexe+i)
+            end do
+
+            ! Local nodal coordinates
+            do i = 1, 8
+               do j = 1, 3
+                  xl(j,i)=co(j,konl(i))
+               end do
+            end do
+
+            ! Loop through gauss points of each element
+            do gp_id = 1, 8
+               ! gauss3d2: hex, 2-point integration (8 integration points)
+               xi = gauss3d2(1,gp_id)
+               et = gauss3d2(2,gp_id)
+               ze = gauss3d2(3,gp_id)
+
+               ! Get the shape function
+               call shape8h(xi,et,ze,xl,xsj,shp,iflag)
+
+               ! Calculate the Gauss point coordinates
+               do k = 1, 3
+                  idx = (el_id - 1)*24 + (gp_id - 1)*3 + k
+                  do l = 1, 8
+                     elem_gp_coord(idx)=elem_gp_coord(idx)
+     &                     +xl(k,l)*shp(4,l)
+                  end do
+               end do
+
+            end do
+
+         end do
+
+         return
+
+      end subroutine getc3d8elementgausspointcoords