From 4d99ec7ebd77973594134da73ddf28aefb3c98c2 Mon Sep 17 00:00:00 2001
From: Ishaan Desai <ishaandesai@gmail.com>
Date: Mon, 29 Sep 2025 17:05:21 +0200
Subject: [PATCH 01/10] Add functionality for handling element meshes, for
 volumetric coupling

---
 adapter/CCXHelpers.c       | 19 ++++++------
 adapter/CCXHelpers.h       | 10 +++++++
 adapter/ConfigReader.cpp   |  4 +++
 adapter/ConfigReader.h     |  1 +
 adapter/PreciceInterface.c | 60 ++++++++++++++++++++++++++++++++++++--
 adapter/PreciceInterface.h | 13 ++++++++-
 6 files changed, 95 insertions(+), 12 deletions(-)

diff --git a/adapter/CCXHelpers.c b/adapter/CCXHelpers.c
index 62dfef43..06f0871f 100644
--- a/adapter/CCXHelpers.c
+++ b/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...
diff --git a/adapter/CCXHelpers.h b/adapter/CCXHelpers.h
index d63aa0a3..e4d0c69e 100644
--- a/adapter/CCXHelpers.h
+++ b/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
diff --git a/adapter/ConfigReader.cpp b/adapter/ConfigReader.cpp
index 66e8d589..06bc4545 100644
--- a/adapter/ConfigReader.cpp
+++ b/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());
     }
diff --git a/adapter/ConfigReader.h b/adapter/ConfigReader.h
index 4044472d..4ffa85b2 100644
--- a/adapter/ConfigReader.h
+++ b/adapter/ConfigReader.h
@@ -13,6 +13,7 @@
 typedef struct InterfaceConfig {
   char * facesMeshName;
   char * nodesMeshName;
+  char * elementsMeshName;
   char * patchName;
   int    map;
   int    numWriteData;
diff --git a/adapter/PreciceInterface.c b/adapter/PreciceInterface.c
index cc3196d9..a48778fd 100644
--- a/adapter/PreciceInterface.c
+++ b/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);
diff --git a/adapter/PreciceInterface.h b/adapter/PreciceInterface.h
index bc343de2..e98ff7e1 100644
--- a/adapter/PreciceInterface.h
+++ b/adapter/PreciceInterface.h
@@ -33,10 +33,14 @@ typedef struct PreciceInterface {
   int          nodeSetID;
   int *        preciceNodeIDs;
   char *       nodesMeshName;
+  char *       elementsMeshName;
 
-  // Interface face elements
+  // Interface elements
   int     numElements;
   int *   elementIDs;
+  int *   elemIPID;
+  double *elemIPCoordinates;
+
   int *   faceIDs;
   double *faceCenterCoordinates;
   int     faceSetID;
@@ -277,6 +281,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

From 02df2f54d5c8376bcb6d1f70ec9eeb8ac8c16def Mon Sep 17 00:00:00 2001
From: Ishaan Desai <ishaandesai@gmail.com>
Date: Mon, 29 Sep 2025 17:07:31 +0200
Subject: [PATCH 02/10] Add functionality to detect quadrature points of meshes
 having either C3D4 or C3D8 elements

---
 getc3d4elementgausspointcoords.f | 86 ++++++++++++++++++++++++++++++
 getc3d8elementgausspointcoords.f | 91 ++++++++++++++++++++++++++++++++
 2 files changed, 177 insertions(+)
 create mode 100755 getc3d4elementgausspointcoords.f
 create mode 100755 getc3d8elementgausspointcoords.f

diff --git a/getc3d4elementgausspointcoords.f b/getc3d4elementgausspointcoords.f
new file mode 100755
index 00000000..bffc43fd
--- /dev/null
+++ b/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
diff --git a/getc3d8elementgausspointcoords.f b/getc3d8elementgausspointcoords.f
new file mode 100755
index 00000000..c36bb2b4
--- /dev/null
+++ b/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

From 67916c22a0012d47e4958a0aa60c8ee0a2b97049 Mon Sep 17 00:00:00 2001
From: Ishaan Desai <ishaandesai@gmail.com>
Date: Mon, 29 Sep 2025 19:17:31 +0200
Subject: [PATCH 03/10] Add the Gauss point extraction functionality files to
 the Makefile, and remove the old file associated with it

---
 Makefile               |   2 +-
 getgausspointscoords.f | 241 -----------------------------------------
 2 files changed, 1 insertion(+), 242 deletions(-)
 delete mode 100644 getgausspointscoords.f

diff --git a/Makefile b/Makefile
index dfab5eca..e146639e 100644
--- a/Makefile
+++ b/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
 
 
 
diff --git a/getgausspointscoords.f b/getgausspointscoords.f
deleted file mode 100644
index 3979491d..00000000
--- a/getgausspointscoords.f
+++ /dev/null
@@ -1,241 +0,0 @@
-!
-!     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.
-!
-      subroutine getgausspointscoords(co,
-     &  iset,istartset,iendset,ipkon,lakon,kon,
-     &  ialset)
-!
-!     calculation and printout of the lift and drag forces
-!
-      implicit none
-!
-      character*8 lakonl,lakon(*)
-!
-      integer konl(20),ifaceq(8,6),nelem,ii,i,j,i1,i2,j1,
-     &  k1,jj,ig,nope,nopes,
-     &  mint2d,ifacet(6,4),ifacew(8,5),iflag,indexe,jface,istartset(*),
-     &  iendset(*),ipkon(*),kon(*),iset,ialset(*),
-     &  fidx
-!
-      real*8 co(3,*),xl(3,20),shp(4,20),xs2(3,7),dvi,f(0:3),
-     &  vkl(0:3,3),t(3,3),div,
-     &  xl2(3,8),xsj2(3),
-     &  shp2(7,8),xi,et,xsj,temp,xi3d,et3d,ze3d,weight,
-     &  xlocal20(3,9,6),xlocal4(3,1,4),xlocal10(3,3,4),xlocal6(3,1,5),
-     &  xlocal15(3,4,5),xlocal8(3,4,6),xlocal8r(3,1,6),pres,
-     &  tf(0:3),tn,tt,dd,coords(3)
-!
-      include "gauss.f"
-      include "xlocal.f"
-!
-      data ifaceq /4,3,2,1,11,10,9,12,
-     &            5,6,7,8,13,14,15,16,
-     &            1,2,6,5,9,18,13,17,
-     &            2,3,7,6,10,19,14,18,
-     &            3,4,8,7,11,20,15,19,
-     &            4,1,5,8,12,17,16,20/
-      data ifacet /1,3,2,7,6,5,
-     &             1,2,4,5,9,8,
-     &             2,3,4,6,10,9,
-     &             1,4,3,8,10,7/
-      data ifacew /1,3,2,9,8,7,0,0,
-     &             4,5,6,10,11,12,0,0,
-     &             1,2,5,4,7,14,10,13,
-     &             2,3,6,5,8,15,11,14,
-     &             4,6,3,1,12,15,9,13/
-      data iflag /3/
-!
-!
-!
-!           initialisierung of the flux
-!     
-            f(0)=0.d0
-			
-!			index for the output array
-			fidx=1
-			
-!     
-!
-            do jj=istartset(iset),iendset(iset)
-!     
-               jface=ialset(jj)
-!     
-               nelem=int(jface/10.d0)
-               ig=jface-10*nelem
-               lakonl=lakon(nelem)
-               indexe=ipkon(nelem)
-!     
-               if(lakonl(4:4).eq.'2') then
-                  nope=20
-                  nopes=8
-               elseif(lakonl(4:4).eq.'8') then
-                  nope=8
-                  nopes=4
-               elseif(lakonl(4:5).eq.'10') then
-                  nope=10
-                  nopes=6
-               elseif(lakonl(4:4).eq.'4') then
-                  nope=4
-                  nopes=3
-               elseif(lakonl(4:5).eq.'15') then
-                  nope=15
-               elseif(lakonl(4:4).eq.'6') then
-                  nope=6
-               endif
-!     
-               if(lakonl(4:5).eq.'8R') then
-                  mint2d=1
-               elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) 
-     &             then
-                  if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'S').or.
-     &                 (lakonl(7:7).eq.'E')) then
-                     mint2d=2
-                  else
-                     mint2d=4
-                  endif
-               elseif(lakonl(4:4).eq.'2') then
-                  mint2d=9
-               elseif(lakonl(4:5).eq.'10') then
-                  mint2d=3
-               elseif(lakonl(4:4).eq.'4') then
-                  mint2d=1
-               endif
-!     
-!     local topology
-!     
-               do i=1,nope
-                  konl(i)=kon(indexe+i)
-               enddo
-!     
-!     computation of the coordinates of the local nodes
-!     
-               do i=1,nope
-                  do j=1,3
-                     xl(j,i)=co(j,konl(i))
-                  enddo
-               enddo
-!     
-!     treatment of wedge faces
-!     
-               if(lakonl(4:4).eq.'6') then
-                  mint2d=1
-                  if(ig.le.2) then
-                     nopes=3
-                  else
-                     nopes=4
-                  endif
-               endif
-               if(lakonl(4:5).eq.'15') then
-                  if(ig.le.2) then
-                     mint2d=3
-                     nopes=6
-                  else
-                     mint2d=4
-                     nopes=8
-                  endif
-               endif
-!     
-               if((nope.eq.20).or.(nope.eq.8)) then
-                  do i=1,nopes
-                     do j=1,3
-                        xl2(j,i)=co(j,konl(ifaceq(i,ig)))
-                     enddo
-                  enddo
-               elseif((nope.eq.10).or.(nope.eq.4)) then
-                  do i=1,nopes
-                     do j=1,3
-                        xl2(j,i)=co(j,konl(ifacet(i,ig)))
-                     enddo
-                  enddo
-               else
-                  do i=1,nopes
-                     do j=1,3
-                        xl2(j,i)=co(j,konl(ifacew(i,ig)))
-                     enddo
-                  enddo
-               endif
-!     
-
-
-               do i=1,mint2d
-!     
-!     local coordinates of the surface integration
-!     point within the surface local coordinate system
-!     
-                  if((lakonl(4:5).eq.'8R').or.
-     &                 ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
-                     xi=gauss2d1(1,i)
-                     et=gauss2d1(2,i)
-                     weight=weight2d1(i)
-                  elseif((lakonl(4:4).eq.'8').or.
-     &                    (lakonl(4:6).eq.'20R').or.
-     &                    ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
-                     xi=gauss2d2(1,i)
-                     et=gauss2d2(2,i)
-                     weight=weight2d2(i)
-                  elseif(lakonl(4:4).eq.'2') then
-                     xi=gauss2d3(1,i)
-                     et=gauss2d3(2,i)
-                     weight=weight2d3(i)
-                  elseif((lakonl(4:5).eq.'10').or.
-     &                    ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
-                     xi=gauss2d5(1,i)
-                     et=gauss2d5(2,i)
-                     weight=weight2d5(i)
-                  elseif((lakonl(4:4).eq.'4').or.
-     &                    ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
-                     xi=gauss2d4(1,i)
-                     et=gauss2d4(2,i)
-                     weight=weight2d4(i)
-                  endif
-!     
-!     local surface normal
-!     
-                  if(nopes.eq.8) then
-                     call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
-                  elseif(nopes.eq.4) then
-                     call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
-                  elseif(nopes.eq.6) then
-                     call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
-                  else
-                     call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
-                  endif
-!
-!                 global coordinates of the integration point
-!
-                  do j1=1,3
-                     coords(j1)=0.d0
-                     do i1=1,nopes
-                        coords(j1)=coords(j1)+shp2(4,i1)*xl2(j1,i1)
-                     enddo
-                  enddo
-                  
-                  print *, coords(1), coords(2), coords(3)
-
-!     
-               enddo
-			   
-			   
-            enddo
-!
-!     
-      return
-      end
-      
-      
-

From d7efaf2a01a257802e9d17dd23b30eff156658e3 Mon Sep 17 00:00:00 2001
From: Ishaan Desai <ishaandesai@gmail.com>
Date: Mon, 29 Sep 2025 19:43:19 +0200
Subject: [PATCH 04/10] Add definition of numIPTotal

---
 adapter/PreciceInterface.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/adapter/PreciceInterface.h b/adapter/PreciceInterface.h
index e98ff7e1..c7477d81 100644
--- a/adapter/PreciceInterface.h
+++ b/adapter/PreciceInterface.h
@@ -35,11 +35,11 @@ typedef struct PreciceInterface {
   char *       nodesMeshName;
   char *       elementsMeshName;
 
-  // Interface elements
   int     numElements;
   int *   elementIDs;
   int *   elemIPID;
   double *elemIPCoordinates;
+  int     numIPTotal;
 
   int *   faceIDs;
   double *faceCenterCoordinates;

From 3893cb5d4f14c26e237b2f468d22271d234ccbc7 Mon Sep 17 00:00:00 2001
From: Ishaan Desai <ishaandesai@gmail.com>
Date: Mon, 29 Sep 2025 19:48:44 +0200
Subject: [PATCH 05/10] Add definition of elementSetID

---
 adapter/PreciceInterface.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/adapter/PreciceInterface.h b/adapter/PreciceInterface.h
index c7477d81..b2aba466 100644
--- a/adapter/PreciceInterface.h
+++ b/adapter/PreciceInterface.h
@@ -38,6 +38,7 @@ typedef struct PreciceInterface {
   int     numElements;
   int *   elementIDs;
   int *   elemIPID;
+  int     elementSetID;
   double *elemIPCoordinates;
   int     numIPTotal;
 

From 96b3f0d1181b5983f5f64d637ce5d5b69c4b82b1 Mon Sep 17 00:00:00 2001
From: Gerasimos Chourdakis <gerasimos.chourdakis@ipvs.uni-stuttgart.de>
Date: Sun, 24 May 2026 19:07:49 +0200
Subject: [PATCH 06/10] Port documentation from the website repository

Original PR: https://github.com/precice/precice.github.io/pull/530

co-authored-by Ishaan Desai <ishaandesai@gmail.com>
---
 docs/adapter-calculix-configure.md | 29 ++++++++++++++++++++++++++++-
 docs/adapter-calculix-overview.md  |  2 +-
 2 files changed, 29 insertions(+), 2 deletions(-)

diff --git a/docs/adapter-calculix-configure.md b/docs/adapter-calculix-configure.md
index 486cd236..8acbefa5 100644
--- a/docs/adapter-calculix-configure.md
+++ b/docs/adapter-calculix-configure.md
@@ -25,11 +25,19 @@ precice-config-file: ../precice-config.xml
 The adapter allows us to use several participants in one simulation (e.g., several instances of Calculix if several solid objects are taken into account). The name of the participant `Calculix` must match the specification of the participant on the command line when running the executable of `CCX` with the adapter being used (this is described later). Also, the name must be the same as the one used in the preCICE configuration file `precice-config.xml`.  
 
 One participant may have several coupling interfaces. Note that each interface specification starts with a dash.
-Depending on the data you need to read and write, the interface should define either a `faces-mesh` (or `mesh` as a synonym) where the data points are centers of faces (computed by the adapter) or a mesh made of CalculiX vertices, with the keyword `nodes-mesh`. An interface made of faces should be defined in the CalculiX case using the `*SURFACE` command, whereas meshes with nodes should define these nodes using `*NSET`. Using the wrong family of mesh (e.g. reading forces on faces) throws an error. If you need both kinds of meshes, you should define more than one interface.
+Depending on the data you need to read and write, the interface should define a mesh of the following types
+
+* a `faces-mesh` (or `mesh` as a synonym) where the data points are centers of faces (computed by the adapter). An interface made of faces should be defined in the CalculiX case using the `*SURFACE` command.
+* a `nodes-mesh` where the data points are the nodal vertices. An interface made of nodes should define these nodes using `*NSET`.
+* a `elements-mesh` where the data points are the quadrature points of the elements of a mesh. The mesh should be defined by nodes using `*NEST`.
+
+Using the wrong family of mesh (e.g. reading forces on faces) throws an error. If you need both kinds of meshes, you should define more than one interface.
 
 For FSI simulations the mesh type of an interface is always `nodes-mesh`, as forces and displacement are defined on nodes. The name of this mesh, `Calculix_Mesh`, must match the mesh name given in the preCICE configuration file. In CHT simulations, `faces-meshes` are usually chosen, as they are needed to apply heat fluxes or convective heat transfer.
 For defining which nodes of the CalculiX domain belong to the FSI interface, a node set needs to be defined in the CalculiX input files. The name of this node set must match the name of the patch (here: "interface").  
 
+For multiscale mechanics simulations, the mesh type is always `elements-mesh`. The stresses, strains, and the material stiffness are defined on the quadrature points.
+
 In the current FSI example, the adapter reads forces from preCICE and feeds displacement deltas (not absolute displacements, but the change of the displacements relative to the last time step) to preCICE. This is defined with the keywords `read-data` and `write-data`, respectively. The names (here: `Forces` and `DisplacementDeltas`) again need to match the specifications in the preCICE configuration file. In the current example, the coupled fluid solver expects displacement deltas instead of displacements. However, the adapter is capable of writing either type. Just use `write-data: [Displacements]` for absolute displacements rather than relative changes being transferred in each time step. Valid `readData` keywords in CalculiX are:
 
 On faces-mesh:
@@ -45,6 +53,18 @@ On nodes-mesh:
 * Displacements (Use `*BOUNDARY`)
 * Temperature (Use `*BOUNDARY`)
 
+On elements-mesh:
+
+* stresses1to3 (components (1,1), (2,2), (3,3) of the stress tensor)
+* stresses4to6 (components (2,3), (1,3), (1,2) of the stress tensor)
+* cmat1 (components (1,1), (1,2), (1,3) of the material stiffness tensor)
+* cmat2 (components (1,4), (1,5), (1,6) of the material stiffness tensor)
+* cmat3 (components (2,2), (2,3), (2,4) of the material stiffness tensor)
+* cmat4 (components (2,5), (2,5), (3,3) of the material stiffness tensor)
+* cmat5 (components (3,4), (3,5), (3,6) of the material stiffness tensor)
+* cmat6 (components (4,4), (4,5), (4,6) of the material stiffness tensor)
+* cmat7 (components (5,5), (5,6), (6,6) of the material stiffness tensor)
+
 Have a look at the CalculiX documentation for a detailed description of each of these commands. There is an [online (but outdated) version](https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node1.html) and an [up-to-date PDF version](http://www.dhondt.de/ccx_2.19.pdf).
 
  Valid `writeData` keywords are:
@@ -70,6 +90,11 @@ From CalculiX version 2.15, additional `writeData` keywords are available:
 * Velocities
 ```
 
+On elements-mesh:
+
+* strains1to3 (components (1,1), (2,2), (3,3) of the strain tensor)
+* strains4to6 (components (2,3), (1,3), (1,2) of the strain tensor)
+
 Note that the square brackets imply that several read- and write-data types can be used on a single interface. This is not needed for FSI simulations (but for CHT simulations). Lastly, the `precice-config-file` needs to be identified including its location. In this example, the file is called `precice-config.xml` and is located one directory above the folder, in which the YAML configuration file lies.
 
 ## CalculiX case input file
@@ -144,6 +169,8 @@ The preCICE CalculiX adapter should support most elements when using `nodes-mesh
 
 When using face meshes, only tetrahedra and hexahedra are supported.
 
+When using `elements-mesh`, linear tetrahedral (C3D4) and hexahedral (C3D8) elements are supported.
+
 ### Coupling to 2D simulations
 
 The adapter supports quasi 2D simulations when the z-direction is ignored. If you set the preCICE interface dimension to 2, the adapter will map data from the CalculiX 3D simulation to 2D space and vice-versa. The 3D simulation should be made of solid elements (or shells) of unit thickness.
diff --git a/docs/adapter-calculix-overview.md b/docs/adapter-calculix-overview.md
index a1f8a60b..f6d2f3ef 100644
--- a/docs/adapter-calculix-overview.md
+++ b/docs/adapter-calculix-overview.md
@@ -21,7 +21,7 @@ The latest supported CalculiX version is {{site.calculix_version}}. If you alrea
 
 ## History
 
-The adapter was initially developed for conjugate heat transfer (CHT) simulations via preCICE by Lucia Cheung in the scope of her master’s thesis[^1], in cooperation with [SimScale](https://www.simscale.com/). For running the adapter for CHT simulations refer to this thesis. The adapter was extended to fluid-structure interaction by Alexander Rusch[^2].
+The adapter was initially developed for conjugate heat transfer (CHT) simulations via preCICE by Lucia Cheung in the scope of her master’s thesis[^1], in cooperation with [SimScale](https://www.simscale.com/). For running the adapter for CHT simulations refer to this thesis. The adapter was extended to fluid-structure interaction by Alexander Rusch[^2]. The adapter was further extended for volumetric coupling of multiscale mechanics problems by Ibrahim Kaleel and Ishaan Desai.
 
 ## References
 

From 060a9b40dfc38f39a0aef323dbdc69984eb21487 Mon Sep 17 00:00:00 2001
From: Ishaan Desai <ishaandesai@gmail.com>
Date: Sat, 1 Aug 2026 17:33:32 +0200
Subject: [PATCH 07/10] Merge the functionality to extract Gauss points for
 C3D4 and C3D8 elements

---
 Makefile                                      |  2 +-
 adapter/PreciceInterface.c                    | 24 +++---
 getc3d4elementgausspointcoords.f              | 86 -------------------
 ...ntcoords.f => getelementgausspointcoords.f | 60 +++++++++----
 4 files changed, 54 insertions(+), 118 deletions(-)
 delete mode 100755 getc3d4elementgausspointcoords.f
 rename getc3d8elementgausspointcoords.f => getelementgausspointcoords.f (55%)
 mode change 100755 => 100644

diff --git a/Makefile b/Makefile
index 331840d0..7800c94b 100644
--- a/Makefile
+++ b/Makefile
@@ -78,7 +78,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 getc3d8elementgausspointcoords.f getc3d4elementgausspointcoords.f
+SCCXF += getflux.f getkdeltatemp.f getelementgausspointcoords.f
 
 
 
diff --git a/adapter/PreciceInterface.c b/adapter/PreciceInterface.c
index f614bac3..e57dc059 100644
--- a/adapter/PreciceInterface.c
+++ b/adapter/PreciceInterface.c
@@ -604,24 +604,24 @@ void PreciceInterface_ConfigureElementsMesh(PreciceInterface *interface, Simulat
 
   enum ElemType elemType = findSimulationMeshType(sim);
 
+  // Gauss point extraction is supported only for tetrahedra and hexahedra elements.
+  int nodesPerElement;
   if (elemType == TETRAHEDRA) {
-    FORTRAN(getc3d4elementgausspointcoords, (&numElements,
-                                             interface->elementIDs,
-                                             sim->co,
-                                             sim->kon,
-                                             sim->ipkon,
-                                             interface->elemIPCoordinates));
+    nodesPerElement = 4;
   } else if (elemType == HEXAHEDRA) {
-    FORTRAN(getc3d8elementgausspointcoords, (&numElements,
-                                             interface->elementIDs,
-                                             sim->co,
-                                             sim->kon,
-                                             sim->ipkon,
-                                             interface->elemIPCoordinates));
+    nodesPerElement = 8;
   } else {
     supportedElementError();
   }
 
+  FORTRAN(getelementgausspointcoords, (&numElements,
+                                       interface->elementIDs,
+                                       &nodesPerElement,
+                                       sim->co,
+                                       sim->kon,
+                                       sim->ipkon,
+                                       interface->elemIPCoordinates));
+
   precicec_setMeshVertices(interface->elementsMeshName, interface->numIPTotal, interface->elemIPCoordinates, interface->elemIPID);
 }
 
diff --git a/getc3d4elementgausspointcoords.f b/getc3d4elementgausspointcoords.f
deleted file mode 100755
index bffc43fd..00000000
--- a/getc3d4elementgausspointcoords.f
+++ /dev/null
@@ -1,86 +0,0 @@
-!
-!     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
diff --git a/getc3d8elementgausspointcoords.f b/getelementgausspointcoords.f
old mode 100755
new mode 100644
similarity index 55%
rename from getc3d8elementgausspointcoords.f
rename to getelementgausspointcoords.f
index c36bb2b4..ccbe107b
--- a/getc3d8elementgausspointcoords.f
+++ b/getelementgausspointcoords.f
@@ -18,7 +18,12 @@
 !
 !     This function is a copy of the function for calculation and printout of the lift and drag forces
 !
-      subroutine getc3d8elementgausspointcoords(nelem, elem_ids,
+!     Computes the physical coordinates of the Gauss/integration points
+!     of a set of elements. Supports either C3D4 (tetrahedra, 1 GP) or
+!     C3D8 (hexahedra, 8 GP) elements, selected via the "nope" argument
+!     (number of nodes per element: 4 or 8).
+!
+      subroutine getelementgausspointcoords(nelem, elem_ids, nope,
      &    co, kon, ipkon, elem_gp_coord)
 
          implicit none
@@ -26,6 +31,7 @@ subroutine getc3d8elementgausspointcoords(nelem, elem_ids,
          !Input/Output variables
          integer     :: nelem            ! Number of elements
          integer     :: elem_ids(*)      ! Element IDs, need to be adjusted for Fortran
+         integer     :: nope             ! Nodes per element: 4 (C3D4) or 8 (C3D8)
          real(8)     :: co(3,*)          ! Nodal coordinates of all nodes
          integer     :: kon(*)
          integer     :: ipkon(*)
@@ -34,15 +40,22 @@ subroutine getc3d8elementgausspointcoords(nelem, elem_ids,
 
          ! Internal variables
          integer(4)     :: iel, el_id, indexe, i, j, k, l
-         integer(4)     :: iflag, idx, konl(8), gp_id
+         integer(4)     :: iflag, idx, konl(8), gp_id, ngp
          real(8)        :: xl(3,8), xi, et, ze, xsj, shp(4,20)
 
          include "gauss.f"
 
          data iflag /1/
 
+         ! Number of Gauss points per element: 1 for tetrahedra, 8 for hexahedra
+         if (nope .eq. 4) then
+            ngp = 1
+         else
+            ngp = 8
+         end if
+
          ! Initialize gauss point coordinates to zero
-         do l = 1, nelem*8*3
+         do l = 1, nelem*ngp*3
             elem_gp_coord(l) = 0.0
          end do
 
@@ -51,32 +64,41 @@ subroutine getc3d8elementgausspointcoords(nelem, elem_ids,
 
             indexe=ipkon(el_id)
 
-            ! connectivity
-            do i = 1, 8
+            ! connectivity + local nodal coordinates
+            do i = 1, nope
                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)
+            do gp_id = 1, ngp
+
+               if (nope .eq. 4) then
+                  ! gauss3d4: tet, 1 integration point
+                  xi = gauss3d4(1,gp_id)
+                  et = gauss3d4(2,gp_id)
+                  ze = gauss3d4(3,gp_id)
+
+                  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
+               else
+                  ! 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)
+                  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
+               end if
 
                ! Calculate the Gauss point coordinates
+               ! Note: indexed by the loop counter (iel), not the element ID
+               ! (el_id), since elem_gp_coord is a compact buffer holding only
+               ! the "nelem" elements of this interface, in the order given by
+               ! elem_ids.
                do k = 1, 3
-                  idx = (el_id - 1)*24 + (gp_id - 1)*3 + k
-                  do l = 1, 8
+                  idx = (iel - 1)*ngp*3 + (gp_id - 1)*3 + k
+                  do l = 1, nope
                      elem_gp_coord(idx)=elem_gp_coord(idx)
      &                     +xl(k,l)*shp(4,l)
                   end do
@@ -88,4 +110,4 @@ subroutine getc3d8elementgausspointcoords(nelem, elem_ids,
 
          return
 
-      end subroutine getc3d8elementgausspointcoords
+      end subroutine getelementgausspointcoords

From c831b7eb423c190d008722683bb3b22d53d287bb Mon Sep 17 00:00:00 2001
From: Ishaan Desai <ishaandesai@gmail.com>
Date: Sat, 1 Aug 2026 17:41:38 +0200
Subject: [PATCH 08/10] Remove incorrectly added parts of the documentation

---
 docs/configure.md | 19 +------------------
 docs/overview.md  |  2 +-
 2 files changed, 2 insertions(+), 19 deletions(-)

diff --git a/docs/configure.md b/docs/configure.md
index 9fae9346..789c0e56 100644
--- a/docs/configure.md
+++ b/docs/configure.md
@@ -31,7 +31,7 @@ Depending on the data you need to read and write, the interface should define a
 
 * a `faces-mesh` (or `mesh` as a synonym) where the data points are centers of faces (computed by the adapter). An interface made of faces should be defined in the CalculiX case using the `*SURFACE` command.
 * a `nodes-mesh` where the data points are the nodal vertices. An interface made of nodes should define these nodes using `*NSET`.
-* a `elements-mesh` where the data points are the quadrature points of the elements of a mesh. The mesh should be defined by nodes using `*NEST`.
+* a `elements-mesh` where the data points are the quadrature points of the elements of a mesh. The mesh should be defined by nodes using `*NEST`. **Note**: `elements-mesh` is still experimental.
 
 Using the wrong family of mesh (e.g. reading forces on faces) throws an error. If you need both kinds of meshes, you should define more than one interface.
 
@@ -55,18 +55,6 @@ On nodes-mesh:
 * Displacements (Use `*BOUNDARY`)
 * Temperature (Use `*BOUNDARY`)
 
-On elements-mesh:
-
-* stresses1to3 (components (1,1), (2,2), (3,3) of the stress tensor)
-* stresses4to6 (components (2,3), (1,3), (1,2) of the stress tensor)
-* cmat1 (components (1,1), (1,2), (1,3) of the material stiffness tensor)
-* cmat2 (components (1,4), (1,5), (1,6) of the material stiffness tensor)
-* cmat3 (components (2,2), (2,3), (2,4) of the material stiffness tensor)
-* cmat4 (components (2,5), (2,5), (3,3) of the material stiffness tensor)
-* cmat5 (components (3,4), (3,5), (3,6) of the material stiffness tensor)
-* cmat6 (components (4,4), (4,5), (4,6) of the material stiffness tensor)
-* cmat7 (components (5,5), (5,6), (6,6) of the material stiffness tensor)
-
 Have a look at the CalculiX documentation for a detailed description of each of these commands. There is an [online (but outdated) version](https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node1.html) and an [up-to-date PDF version](http://www.dhondt.de/ccx_2.19.pdf).
 
  Valid `writeData` keywords are:
@@ -92,11 +80,6 @@ From CalculiX version 2.15, additional `writeData` keywords are available:
 * Velocities
 ```
 
-On elements-mesh:
-
-* strains1to3 (components (1,1), (2,2), (3,3) of the strain tensor)
-* strains4to6 (components (2,3), (1,3), (1,2) of the strain tensor)
-
 Note that the square brackets imply that several read- and write-data types can be used on a single interface. This is not needed for FSI simulations (but for CHT simulations). Lastly, the `precice-config-file` needs to be identified including its location. In this example, the file is called `precice-config.xml` and is located one directory above the folder, in which the YAML configuration file lies.
 
 ## CalculiX case input file
diff --git a/docs/overview.md b/docs/overview.md
index 5c11e103..08422a6e 100644
--- a/docs/overview.md
+++ b/docs/overview.md
@@ -32,7 +32,7 @@ Compatibility with preCICE:
 
 ## History
 
-The adapter was initially developed for conjugate heat transfer (CHT) simulations via preCICE by Lucia Cheung in the scope of her master’s thesis[^1], in cooperation with [SimScale](https://www.simscale.com/). For running the adapter for CHT simulations refer to this thesis. The adapter was extended to fluid-structure interaction by Alexander Rusch[^2]. The adapter was further extended for volumetric coupling of multiscale mechanics problems by Ibrahim Kaleel and Ishaan Desai.
+The adapter was initially developed for conjugate heat transfer (CHT) simulations via preCICE by Lucia Cheung in the scope of her master’s thesis[^1], in cooperation with [SimScale](https://www.simscale.com/). For running the adapter for CHT simulations refer to this thesis. The adapter was extended to fluid-structure interaction by Alexander Rusch[^2]. The adapter was further extended for volumetric coupling by Ibrahim Kaleel and Ishaan Desai.
 
 ## References
 

From 857d5c80ff87988ae081c9ec3c69629e3c7df31c Mon Sep 17 00:00:00 2001
From: Ishaan Desai <ishaandesai@gmail.com>
Date: Sat, 1 Aug 2026 17:49:18 +0200
Subject: [PATCH 09/10] Add warning for elements mesh support

---
 adapter/PreciceInterface.c | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/adapter/PreciceInterface.c b/adapter/PreciceInterface.c
index e57dc059..fd324caf 100644
--- a/adapter/PreciceInterface.c
+++ b/adapter/PreciceInterface.c
@@ -585,6 +585,9 @@ static enum ElemType findSimulationMeshType(SimulationData *sim)
 
 void PreciceInterface_ConfigureElementsMesh(PreciceInterface *interface, SimulationData *sim)
 {
+  printf("WARNING: Elements-mesh support is experimental. Use with caution.\n");
+  fflush(stdout);
+
   char *elementSetName    = interface->name;
   interface->elementSetID = getSetID(elementSetName, sim->set, sim->nset);
   interface->numElements  = getNumSetElements(interface->elementSetID, sim->istartset, sim->iendset);

From 7b331cbc9a4792c21177dc76f196334d111a887d Mon Sep 17 00:00:00 2001
From: Ishaan Desai <ishaandesai@gmail.com>
Date: Sat, 1 Aug 2026 22:21:58 +0200
Subject: [PATCH 10/10] Add CHANGELOG entry

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index ed7fcee4..b48fd2e1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/).
 
 ### Added
 
+- Added preliminary functionality to support volumetric coupling by extracting integration points from an element mesh ([#146](https://github.com/precice/calculix-adapter/pull/146)).
 - Extended velocities writing for the 2D3D case ([#137](https://github.com/precice/calculix-adapter/pull/137)).
 - Added Ubuntu 24.04 to the CI ([#140](https://github.com/precice/calculix-adapter/pull/140)).