diff --git a/Makefile b/Makefile index a201b03e..b1bf7b6d 100644 --- a/Makefile +++ b/Makefile @@ -62,7 +62,7 @@ else CC = mpicc endif -FFLAGS = -Wall -O3 -fopenmp $(INCLUDES) ${ADDITIONAL_FFLAGS} +FFLAGS = -Wall -O3 -fopenmp -fallow-argument-mismatch $(INCLUDES) ${ADDITIONAL_FFLAGS} # Note for GCC 10 or newer: add -fallow-argument-mismatch in the above flags FC = mpifort # FC = mpif90 @@ -73,8 +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 getc3d4elementgausspointcoords.f diff --git a/adapter/CCXHelpers.c b/adapter/CCXHelpers.c index 62dfef43..34858a88 100644 --- a/adapter/CCXHelpers.c +++ b/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) { @@ -409,6 +420,18 @@ int getXloadIndexOffset(enum xloadVariable xloadVar) } } +void getElementStrain(int strainIdx, int numIPTotal, double *eei, double *strainData) +{ + int i, idx; + // Loop through all element and respective gauss points + for (i = 0; i < numIPTotal; i++) { + idx = i * 6 + strainIdx; //TODO: Add explanation for 6 + strainData[i * 3] = eei[idx]; + strainData[i * 3 + 1] = eei[idx + 1]; + strainData[i * 3 + 2] = eei[idx + 2]; + } +} + void setXload(double *xload, int *xloadIndices, double *values, int numValues, enum xloadVariable xloadVar) { ITG i; @@ -470,6 +493,30 @@ void setNodeDisplacements(double *displacements, ITG numNodes, int dim, int *xbo } } +void setElementStiffness(int stiffnessIdx, int numIPTotal, double *stiffnessData, double *xstiff) +{ + int i, idx; + // Loop through all element and respective gauss points + for (i = 0; i < numIPTotal; i++) { + idx = i * 27 + stiffnessIdx; //TODO: Add explanation for 27 + xstiff[idx] = stiffnessData[i * 3]; + xstiff[idx + 1] = stiffnessData[i * 3 + 1]; + xstiff[idx + 2] = stiffnessData[i * 3 + 2]; + } +} + +void setElementStress(int stressIdx, int numIPTotal, double *stressData, double *stx) +{ + int i, idx; + // Loop through all element and respective gauss points + for (i = 0; i < numIPTotal; i++) { + idx = i * 6 + stressIdx; //TODO: Add explanation for 6 + stx[idx] = stressData[i * 3]; + stx[idx + 1] = stressData[i * 3 + 1]; + stx[idx + 2] = stressData[i * 3 + 2]; + } +} + bool isSteadyStateSimulation(ITG *nmethod) { return *nmethod == 1; diff --git a/adapter/CCXHelpers.h b/adapter/CCXHelpers.h index d63aa0a3..47e3550f 100644 --- a/adapter/CCXHelpers.h +++ b/adapter/CCXHelpers.h @@ -51,7 +51,18 @@ enum CouplingDataType { TEMPERATURE, DISPLACEMENTS, DISPLACEMENTDELTAS, VELOCITIES, - POSITIONS }; + POSITIONS, + STRAIN1TO3, + STRAIN4TO6, + STRESS1TO3, + STRESS4TO6, + CMAT1, + CMAT2, + CMAT3, + CMAT4, + CMAT5, + CMAT6, + CMAT7 }; /** * @brief Type of element used for faces mesh, where we assume only one type of element is used. @@ -102,6 +113,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 +262,15 @@ 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 numIPTotal: 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 numIPTotal, double *eei, double *strainData); + /** * @brief Modifies the values of a DFLUX or FILM boundary condition * @param xload: CalculiX array for the loads @@ -317,6 +347,18 @@ void setNodeForces(double *forces, ITG numNodes, int dim, int *xforcIndices, dou */ void setNodeDisplacements(double *displacements, ITG numNodes, int dim, int *xbounIndices, double *xboun); +/** + * @brief Modifies the values of stress of the elements at the interface + * @param TODO + */ +void setElementStiffness(int stiffnessIdx, int numIPTotal, double *stiffnessData, double *xstiff); + +/** + * @brief Modifies the values of stress of the elements at the interface + * @param TODO + */ +void setElementStress(int stressIdx, int numIPTotal, double *stressData, double *stx); + /** * @brief Returns whether it is a steady-state simulation based on the value of nmethod * @param nmethod: CalculiX variable with information regarding the type of analysis diff --git a/adapter/ConfigReader.cpp b/adapter/ConfigReader.cpp index e49bc5df..7c4f857f 100644 --- a/adapter/ConfigReader.cpp +++ b/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().c_str()); } + if (config["participants"][participantName]["interfaces"][i]["elements"]) { + interface.elementsMeshName = strdup(config["participants"][participantName]["interfaces"][i]["elements"].as().c_str()); + } + std::string patchName = config["participants"][participantName]["interfaces"][i]["patch"].as(); std::transform(patchName.begin(), patchName.end(), patchName.begin(), toupper); interface.patchName = strdup(patchName.c_str()); diff --git a/adapter/ConfigReader.h b/adapter/ConfigReader.h old mode 100644 new mode 100755 index 4044472d..4ffa85b2 --- 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..1590b057 100644 --- a/adapter/PreciceInterface.c +++ b/adapter/PreciceInterface.c @@ -4,6 +4,7 @@ * Heat transfer adapter developed by Lucía Cheung with the support of SimScale GmbH * * * * Adapter extended to fluid-structure interaction by Alexander Rusch * + * Adapter extended to multiscale mechanics by Ibrahim Kaleel and Ishaan Desai * * * *********************************************************************************************/ @@ -61,7 +62,7 @@ void Precice_Setup(char *configFilename, char *participantName, SimulationData * // Initialize coupling data printf("Initializing coupling data\n"); fflush(stdout); - Precice_ReadCouplingData(sim); + //Precice_ReadCouplingData(sim); } void Precice_AdjustSolverTimestep(SimulationData *sim) @@ -73,9 +74,9 @@ void Precice_AdjustSolverTimestep(SimulationData *sim) fflush(stdout); // For steady-state simulations, we will always compute the converged steady-state solution in one coupling step - *sim->theta = 0; - *sim->tper = 1; - *sim->dtheta = 1; + //*sim->theta = 0; + //*sim->tper = 1; + //*sim->dtheta = 1; // Set the solver time step to be the same as the coupling time step sim->solver_dt = precice_dt; @@ -264,6 +265,74 @@ void Precice_ReadCouplingData(SimulationData *sim) } printf("Reading DISPLACEMENTS coupling data.\n"); break; + case CMAT1: + // READ MATERIAL MATRIX COMPONENTS - C11, C12, C13 + precicec_readData(interfaces[i]->couplingMeshName, interfaces[i]->materialTangent1Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, sim->solver_dt, interfaces[i]->elementIPVectorData); + setElementStiffness(0, interfaces[i]->numIPTotal, interfaces[i]->elementIPVectorData, sim->xstiff); + printf("Reading MATERIAL TANGENT 1 coupling data.\n"); + break; + case CMAT2: + // READ MATERIAL MATRIX COMPONENTS - C14, C15, C16 + precicec_readData(interfaces[i]->couplingMeshName, interfaces[i]->materialTangent2Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, sim->solver_dt, interfaces[i]->elementIPVectorData); + setElementStiffness(3, interfaces[i]->numIPTotal, interfaces[i]->elementIPVectorData, sim->xstiff); + printf("Reading MATERIAL TANGENT 2 coupling data.\n"); + break; + case CMAT3: + // READ MATERIAL MATRIX COMPONENTS - C22, C23, C24 + precicec_readData(interfaces[i]->couplingMeshName, interfaces[i]->materialTangent3Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, sim->solver_dt, interfaces[i]->elementIPVectorData); + setElementStiffness(6, interfaces[i]->numIPTotal, interfaces[i]->elementIPVectorData, sim->xstiff); + printf("Reading MATERIAL TANGENT 3 coupling data.\n"); + break; + case CMAT4: + // READ MATERIAL MATRIX COMPONENTS - C25, C26, C33 + precicec_readData(interfaces[i]->couplingMeshName, interfaces[i]->materialTangent4Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, sim->solver_dt, interfaces[i]->elementIPVectorData); + setElementStiffness(9, interfaces[i]->numIPTotal, interfaces[i]->elementIPVectorData, sim->xstiff); + printf("Reading MATERIAL TANGENT 4 coupling data.\n"); + break; + case CMAT5: + // READ MATERIAL MATRIX COMPONENTS - C34, C35, C36 + precicec_readData(interfaces[i]->couplingMeshName, interfaces[i]->materialTangent5Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, sim->solver_dt, interfaces[i]->elementIPVectorData); + setElementStiffness(12, interfaces[i]->numIPTotal, interfaces[i]->elementIPVectorData, sim->xstiff); + printf("Reading MATERIAL TANGENT 5 coupling data.\n"); + break; + case CMAT6: + // READ MATERIAL MATRIX COMPONENTS - C44, C45, C46 + precicec_readData(interfaces[i]->couplingMeshName, interfaces[i]->materialTangent6Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, sim->solver_dt, interfaces[i]->elementIPVectorData); + setElementStiffness(15, interfaces[i]->numIPTotal, interfaces[i]->elementIPVectorData, sim->xstiff); + printf("Reading MATERIAL TANGENT 6 coupling data.\n"); + break; + case CMAT7: + // READ MATERIAL MATRIX COMPONENTS - C55, C56, C66 + precicec_readData(interfaces[i]->couplingMeshName, interfaces[i]->materialTangent7Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, sim->solver_dt, interfaces[i]->elementIPVectorData); + for (int k = 0; k < interfaces[i]->numIPTotal; k++) { + printf("CMAT7 read from preCICE: %f, %f, %f\n", interfaces[i]->elementIPVectorData[k * 3], interfaces[i]->elementIPVectorData[k * 3 + 1], interfaces[i]->elementIPVectorData[k * 3 + 2]); + } + setElementStiffness(18, interfaces[i]->numIPTotal, interfaces[i]->elementIPVectorData, sim->xstiff); + printf("Reading MATERIAL TANGENT 7 coupling data.\n"); + break; + case STRESS1TO3: + // READ STRESS COMPONENTS - S11, S22, S33 + precicec_readData(interfaces[i]->couplingMeshName, interfaces[i]->stress1to3Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, sim->solver_dt, interfaces[i]->elementIPVectorData); + // for (int k = 0; k < interfaces[i]->numIPTotal; k++) { + // printf("Stresses1to3 read from preCICE: %f, %f, %f\n", interfaces[i]->elementIPVectorData[k * 3], interfaces[i]->elementIPVectorData[k * 3 + 1], interfaces[i]->elementIPVectorData[k * 3 + 2]); + // } + setElementStress(0, interfaces[i]->numIPTotal, interfaces[i]->elementIPVectorData, sim->stx); + // int idxx; + // for (int k = 0; k < interfaces[i]->numIPTotal; k++) { + // idxx = k * 6; + // printf("stx for %d = [%f,%f,%f]\n", idxx, sim->stx[idxx], sim->stx[idxx + 1], sim->stx[idxx + 2]); + // } + printf("Reading STRESS1TO3 coupling data.\n"); + break; + case STRESS4TO6: + // READ STRESS COMPONENTS - S23, S13, S12 + precicec_readData(interfaces[i]->couplingMeshName, interfaces[i]->stress4to6Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, sim->solver_dt, interfaces[i]->elementIPVectorData); + // for (int k = 0; k < interfaces[i]->numIPTotal; k++) { + // printf("Stresses4to6 read from preCICE: %f, %f, %f\n", interfaces[i]->elementIPVectorData[k * 3], interfaces[i]->elementIPVectorData[k * 3 + 1], interfaces[i]->elementIPVectorData[k * 3 + 2]); + // } + setElementStress(3, interfaces[i]->numIPTotal, interfaces[i]->elementIPVectorData, sim->stx); + printf("Reading STRESS4TO6 coupling data.\n"); + break; case DISPLACEMENTDELTAS: printf("DisplacementDeltas cannot be used as read data\n"); fflush(stdout); @@ -279,6 +348,16 @@ void Precice_ReadCouplingData(SimulationData *sim) fflush(stdout); exit(EXIT_FAILURE); break; + case STRAIN1TO3: + printf("Strain 1to3 cannot be used as read data.\n"); + fflush(stdout); + exit(EXIT_FAILURE); + break; + case STRAIN4TO6: + printf("Strain 4to6 cannot be used as read data.\n"); + fflush(stdout); + exit(EXIT_FAILURE); + break; } } } @@ -315,7 +394,7 @@ void Precice_WriteCouplingData(SimulationData *sim) sim->istartset, sim->iendset, sim->ipkon, - sim->lakon, + *sim->lakon, sim->kon, sim->ialset, sim->ielmat, @@ -413,6 +492,16 @@ void Precice_WriteCouplingData(SimulationData *sim) precicec_writeData(interfaces[i]->couplingMeshName, interfaces[i]->forces, interfaces[i]->numNodes, interfaces[i]->preciceNodeIDs, interfaces[i]->nodeVectorData); printf("Writing FORCES coupling data.\n"); break; + case STRAIN1TO3: + getElementStrain(0, interfaces[i]->numIPTotal, sim->eei, interfaces[i]->elementIPVectorData); + precicec_writeData(interfaces[i]->couplingMeshName, interfaces[i]->strain1to3Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, interfaces[i]->elementIPVectorData); + printf("Writing STRAIN1TO3 coupling data.\n"); + break; + case STRAIN4TO6: + getElementStrain(3, interfaces[i]->numIPTotal, sim->eei, interfaces[i]->elementIPVectorData); + precicec_writeData(interfaces[i]->couplingMeshName, interfaces[i]->strain4to6Data, interfaces[i]->numIPTotal, interfaces[i]->elemIPID, interfaces[i]->elementIPVectorData); + printf("Writing STRAIN4TO6 coupling data.\n"); + break; } } // Cleanup data @@ -450,8 +539,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 mesh, nodes mesh or an element mesh. Check the adapter configuration file.\n"); exit(EXIT_FAILURE); } if (config->nodesMeshName && config->facesMeshName) { @@ -491,6 +580,11 @@ void PreciceInterface_Create(PreciceInterface *interface, SimulationData *sim, I interface->writeData = NULL; interface->readData = NULL; + // Initialize element data points as NULL + interface->elemIPID = NULL; + interface->elemIPCoordinates = NULL; + interface->elementIPVectorData = NULL; + // Initialize preCICE mesh name as NULL interface->couplingMeshName = NULL; @@ -507,6 +601,17 @@ void PreciceInterface_Create(PreciceInterface *interface, SimulationData *sim, I interface->velocities = NULL; interface->forces = NULL; interface->pressure = NULL; + interface->strain1to3Data = NULL; + interface->strain4to6Data = NULL; + interface->stress1to3Data = NULL; + interface->stress4to6Data = NULL; + interface->materialTangent1Data = NULL; + interface->materialTangent2Data = NULL; + interface->materialTangent3Data = NULL; + interface->materialTangent4Data = NULL; + interface->materialTangent5Data = NULL; + interface->materialTangent6Data = NULL; + interface->materialTangent7Data = NULL; // Check if quasi 2D-3D coupling needs to be implemented if (interface->dim == 2) { @@ -545,12 +650,20 @@ 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); } static enum ElemType findSimulationMeshType(SimulationData *sim) { - // Assuming only tetrahedra are used, or only hexaedral, for faces meshes. + // Assuming only tetrahedra or hexahedra are used. // Return first non-zero mesh type. // lakon tab takes 8 chars per element @@ -563,10 +676,57 @@ static enum ElemType findSimulationMeshType(SimulationData *sim) } lakon_ptr += 8; } - return INVALID_ELEMENT; } +void PreciceInterface_ConfigureElementsMesh(PreciceInterface *interface, SimulationData *sim) +{ + printf("Entering ConfigureElementsMesh \n"); + 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(); + } + + // debugging TODO: remove + for (int i = 0; i < interface->numIPTotal; i++) { + printf("Gauss point coordinates: %f, %f, %f\n", interface->elemIPCoordinates[3 * i], interface->elemIPCoordinates[3 * i + 1], interface->elemIPCoordinates[3 * i + 2]); + } + + precicec_setMeshVertices(interface->elementsMeshName, interface->numIPTotal, interface->elemIPCoordinates, interface->elemIPID); +} + void PreciceInterface_ConfigureFaceCentersMesh(PreciceInterface *interface, SimulationData *sim) { // printf("Entering ConfigureFaceCentersMesh \n"); @@ -689,6 +849,8 @@ void PreciceInterface_ConfigureCouplingData(PreciceInterface *interface, Simulat interface->nodeVectorData = malloc(interface->numNodes * 3 * sizeof(double)); interface->faceCenterData = malloc(interface->numElements * sizeof(double)); + interface->elementIPVectorData = malloc(interface->numIPTotal * 3 * sizeof(double)); + // Configure all the read data, then the write data int i; interface->numReadData = config->numReadData; @@ -741,7 +903,52 @@ void PreciceInterface_ConfigureCouplingData(PreciceInterface *interface, Simulat interface->xbounIndices = malloc(interface->numNodes * 3 * sizeof(int)); interface->displacements = strdup(config->readDataNames[i]); getXbounIndices(interface->nodeIDs, interface->numNodes, sim->nboun, sim->ikboun, sim->ilboun, interface->xbounIndices, DISPLACEMENTS); - printf("Read data '%s' found.\n", interface->displacements); + printf("Read data '%s' found.\n", config->readDataNames[i]); + } else if (startsWith(config->readDataNames[i], "cmat1")) { + PreciceInterface_EnsureValidRead(sim, CMAT1); + interface->readData[i] = CMAT1; + interface->materialTangent1Data = strdup(config->readDataNames[i]); + printf("Read data '%s' found.\n", config->readDataNames[i]); + } else if (startsWith(config->readDataNames[i], "cmat2")) { + PreciceInterface_EnsureValidRead(sim, CMAT2); + interface->readData[i] = CMAT2; + interface->materialTangent2Data = strdup(config->readDataNames[i]); + printf("Read data '%s' found.\n", config->readDataNames[i]); + } else if (startsWith(config->readDataNames[i], "cmat3")) { + PreciceInterface_EnsureValidRead(sim, CMAT3); + interface->readData[i] = CMAT3; + interface->materialTangent3Data = strdup(config->readDataNames[i]); + printf("Read data '%s' found.\n", config->readDataNames[i]); + } else if (startsWith(config->readDataNames[i], "cmat4")) { + PreciceInterface_EnsureValidRead(sim, CMAT4); + interface->readData[i] = CMAT4; + interface->materialTangent4Data = strdup(config->readDataNames[i]); + printf("Read data '%s' found.\n", config->readDataNames[i]); + } else if (startsWith(config->readDataNames[i], "cmat5")) { + PreciceInterface_EnsureValidRead(sim, CMAT5); + interface->readData[i] = CMAT5; + interface->materialTangent5Data = strdup(config->readDataNames[i]); + printf("Read data '%s' found.\n", config->readDataNames[i]); + } else if (startsWith(config->readDataNames[i], "cmat6")) { + PreciceInterface_EnsureValidRead(sim, CMAT6); + interface->readData[i] = CMAT6; + interface->materialTangent6Data = strdup(config->readDataNames[i]); + printf("Read data '%s' found.\n", config->readDataNames[i]); + } else if (startsWith(config->readDataNames[i], "cmat7")) { + PreciceInterface_EnsureValidRead(sim, CMAT7); + interface->readData[i] = CMAT7; + interface->materialTangent7Data = strdup(config->readDataNames[i]); + printf("Read data '%s' found.\n", config->readDataNames[i]); + } else if (startsWith(config->readDataNames[i], "stresses1to3")) { + PreciceInterface_EnsureValidRead(sim, STRESS1TO3); + interface->readData[i] = STRESS1TO3; + interface->stress1to3Data = strdup(config->readDataNames[i]); + printf("Read data '%s' found.\n", config->readDataNames[i]); + } else if (startsWith(config->readDataNames[i], "stresses4to6")) { + PreciceInterface_EnsureValidRead(sim, STRESS4TO6); + interface->readData[i] = STRESS4TO6; + interface->stress4to6Data = strdup(config->readDataNames[i]); + printf("Read data '%s' found.\n", config->readDataNames[i]); } else { printf("ERROR: Read data '%s' is not of a known type for the CalculiX-preCICE adapter. Check the adapter configuration file.\n", config->readDataNames[i]); exit(EXIT_FAILURE); @@ -791,6 +998,14 @@ void PreciceInterface_ConfigureCouplingData(PreciceInterface *interface, Simulat interface->writeData[i] = FORCES; interface->forces = strdup(config->writeDataNames[i]); printf("Write data '%s' found.\n", interface->forces); + } else if (isEqual(config->writeDataNames[i], "strains1to3")) { + interface->writeData[i] = STRAIN1TO3; + interface->strain1to3Data = strdup(config->writeDataNames[i]); + printf("Write data '%s' found.\n", config->writeDataNames[i]); + } else if (isEqual(config->writeDataNames[i], "strains4to6")) { + interface->writeData[i] = STRAIN4TO6; + interface->strain4to6Data = strdup(config->writeDataNames[i]); + printf("Write data '%s' found.\n", config->writeDataNames[i]); } else { printf("ERROR: Write data '%s' is not of a known type for the CalculiX-preCICE adapter. Check the adapter configuration file.\n", config->writeDataNames[i]); exit(EXIT_FAILURE); @@ -804,6 +1019,7 @@ void PreciceInterface_FreeData(PreciceInterface *preciceInterface) free(preciceInterface->writeData); free(preciceInterface->elementIDs); free(preciceInterface->faceIDs); + free(preciceInterface->nodeIDs); free(preciceInterface->preciceFaceCenterIDs); free(preciceInterface->faceCenterCoordinates); free(preciceInterface->nodeCoordinates); @@ -819,9 +1035,15 @@ void PreciceInterface_FreeData(PreciceInterface *preciceInterface) freeMapping(preciceInterface->mappingQuasi2D3D); + // Volumteric element related data + free(preciceInterface->elemIPCoordinates); + free(preciceInterface->elemIPID); + free(preciceInterface->elementIPVectorData); + // Mesh names free(preciceInterface->faceCentersMeshName); free(preciceInterface->nodesMeshName); + free(preciceInterface->elementsMeshName); // Data names free(preciceInterface->displacementDeltas); @@ -836,4 +1058,15 @@ void PreciceInterface_FreeData(PreciceInterface *preciceInterface) free(preciceInterface->pressure); free(preciceInterface->temperature); free(preciceInterface->velocities); + free(preciceInterface->strain1to3Data); + free(preciceInterface->strain4to6Data); + free(preciceInterface->stress1to3Data); + free(preciceInterface->stress4to6Data); + free(preciceInterface->materialTangent1Data); + free(preciceInterface->materialTangent2Data); + free(preciceInterface->materialTangent3Data); + free(preciceInterface->materialTangent4Data); + free(preciceInterface->materialTangent5Data); + free(preciceInterface->materialTangent6Data); + free(preciceInterface->materialTangent7Data); } diff --git a/adapter/PreciceInterface.h b/adapter/PreciceInterface.h index bc343de2..872b44b4 100644 --- a/adapter/PreciceInterface.h +++ b/adapter/PreciceInterface.h @@ -43,12 +43,20 @@ typedef struct PreciceInterface { char * faceCentersMeshName; int * preciceFaceCenterIDs; + // Interface volumetric elements + char * elementsMeshName; + int elementSetID; + int numIPTotal; + double *elemIPCoordinates; + int * elemIPID; + // Arrays to store the coupling data double *nodeScalarData; double *node2DScalarData; // Scalar quantities in 2D in case quasi 2D-3D coupling is done double *nodeVectorData; // Forces, displacements, velocities, positions and displacementDeltas are vector quantities double *node2DVectorData; // Vector quantities in 2D in case quasi 2D-3D coupling is done double *faceCenterData; + double *elementIPVectorData; // Vector quantities at the integration points // preCICE mesh name char *couplingMeshName; @@ -66,6 +74,17 @@ typedef struct PreciceInterface { char *velocities; char *forces; char *pressure; + char *strain1to3Data; + char *strain4to6Data; + char *stress1to3Data; + char *stress4to6Data; + char *materialTangent1Data; + char *materialTangent2Data; + char *materialTangent3Data; + char *materialTangent4Data; + char *materialTangent5Data; + char *materialTangent6Data; + char *materialTangent7Data; // Indices that indicate where to apply the boundary conditions / forces int *xloadIndices; @@ -94,6 +113,7 @@ typedef struct PreciceInterface { typedef struct SimulationData { // CalculiX data + ITG * iset; ITG * ialset; ITG * ielmat; ITG * istartset; @@ -130,6 +150,11 @@ typedef struct SimulationData { double *cocon; ITG * ncocon; ITG * mi; + ITG * nea; // element bounds in each thread - start + ITG * neb; // element bounds in each thread - end + double *eei; // Strain tensor values + double *stx; // Stress tensor values + double *xstiff; // Stiffness matrix values // Interfaces int numPreciceInterfaces; @@ -277,6 +302,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 + */ +void PreciceInterface_ConfigureElementsMesh(PreciceInterface *interface, SimulationData *sim); + /** * @brief Terminate execution if the nodes mesh ID is not valid * @param interface diff --git a/ccx_2.20.c b/ccx_2.20.c index 84466232..74a02694 100644 --- a/ccx_2.20.c +++ b/ccx_2.20.c @@ -1308,7 +1308,7 @@ int main(int argc, char *argv[]) /* nmethod=15: Crack propagation */ /* nmethod=16: Feasible direction based on sensitivity information */ if (preciceUsed) { - int isStaticOrDynamic = ((nmethod == 1) || (nmethod == 4)) && (iperturb[0] > 1); + int isStaticOrDynamic = ((nmethod == 1) || (nmethod == 4)) ;//&& (iperturb[0] > 1); Necessary to trigger isStaticorDynamic int isDynamic = ((nmethod == 4) && (iperturb[0] > 1)); int isThermalAnalysis = ithermal[0] >= 2; int isModalDynamic = ((nmethod == 4) && (iperturb[0] < 2)); @@ -1432,6 +1432,55 @@ int main(int argc, char *argv[]) printf("ERROR: This simulation type is not available with preCICE"); exit(0); } + + } else if (isStaticOrDynamic) { + + mpcinfo[0] = memmpc_; + mpcinfo[1] = mpcfree; + mpcinfo[2] = icascade; + mpcinfo[3] = maxlenmpc; + + if (icascade != 0) { + printf(" *ERROR in CalculiX: the matrix structure may"); + printf(" change due to nonlinear equations;"); + printf(" a purely linear calculation is not"); + printf(" feasible; use NLGEOM on the *STEP card."); + FORTRAN(stop, ()); + } + + printf("Starting multiscale linear static analysis via preCICE...\n"); + + linstatic_precice(co, &nk, &kon, &ipkon, &lakon, &ne, nodeboun, ndirboun, xboun, + &nboun, + ipompc, nodempc, coefmpc, labmpc, &nmpc, nodeforc, ndirforc, xforc, + &nforc, nelemload, sideload, xload, &nload, + nactdof, &icol, jq, &irow, neq, &nzl, &nmethod, ikmpc, + ilmpc, ikboun, ilboun, elcon, nelcon, rhcon, nrhcon, + alcon, nalcon, alzero, &ielmat, &ielorien, &norien, orab, &ntmat_, + t0, t1, t1old, ithermal, prestr, &iprestr, vold, iperturb, sti, nzs, + &kode, filab, eme, &iexpl, plicon, + nplicon, plkcon, nplkcon, &xstate, &npmat_, matname, + &isolver, mi, &ncmat_, &nstate_, cs, &mcs, &nkon, &ener, + xbounold, xforcold, xloadold, amname, amta, namta, + &nam, iamforc, iamload, iamt1, iamboun, &ttime, + output, set, &nset, istartset, iendset, ialset, &nprint, prlab, + prset, &nener, trab, inotr, &ntrans, fmpc, ipobody, ibody, xbody, + &nbody, + xbodyold, timepar, thicke, jobnamec, tieset, &ntie, &istep, &nmat, + ielprop, prop, typeboun, &mortar, mpcinfo, tietol, ics, + orname, itempuser, t0g, t1g, + /* PreCICE args */ + preciceParticipantName, configFilename); + + for (i = 0; i < 3; i++) { + nzsprevstep[i] = nzs[i]; + } + + memmpc_ = mpcinfo[0]; + mpcfree = mpcinfo[1]; + icascade = mpcinfo[2]; + maxlenmpc = mpcinfo[3]; + } else if (isStaticNLGEOM) { printf("Starting STATIC analysis via preCICE...\n"); @@ -1506,6 +1555,7 @@ int main(int argc, char *argv[]) t0g, t1g, preciceParticipantName, configFilename); } else { + //printf("DEBUG: nmethod=%d, iperturb[0]=%d, iperturb[1]=%d\n", nmethod, iperturb[0], iperturb[1]); printf("ERROR: Only thermal coupling or FSI is available with preCICE"); exit(0); } 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 diff --git a/linstatic_precice.c b/linstatic_precice.c new file mode 100644 index 00000000..57565041 --- /dev/null +++ b/linstatic_precice.c @@ -0,0 +1,1181 @@ +/* CalculiX - A 3-dimensional finite element program */ +/* Copyright (C) 1998-2023 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. */ + +#include +#include +#include +#include "CalculiX.h" +#ifdef SPOOLES +#include "spooles.h" +#endif +#ifdef SGI +#include "sgi.h" +#endif +#ifdef TAUCS +#include "tau.h" +#endif +#ifdef PARDISO +#include "pardiso.h" +#endif +#ifdef PASTIX +#include "pastix.h" +#endif + +/* Adapter: Add header */ +#include "adapter/PreciceInterface.h" + +void linstatic_precice(double *co, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp, + ITG *ne, + ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, + ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc, + ITG *nmpc, + ITG *nodeforc, ITG *ndirforc, double *xforc, ITG *nforc, + ITG *nelemload, char *sideload, double *xload, + ITG *nload, ITG *nactdof, + ITG **icolp, ITG *jq, ITG **irowp, ITG *neq, ITG *nzl, + ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, + ITG *ilboun, + double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon, + double *alcon, ITG *nalcon, double *alzero, ITG **ielmatp, + ITG **ielorienp, ITG *norien, double *orab, ITG *ntmat_, + double *t0, double *t1, double *t1old, + ITG *ithermal, double *prestr, ITG *iprestr, + double *vold, ITG *iperturb, double *sti, ITG *nzs, + ITG *kode, char *filab, double *eme, + ITG *iexpl, double *plicon, ITG *nplicon, double *plkcon, + ITG *nplkcon, + double **xstatep, ITG *npmat_, char *matname, ITG *isolver, + ITG *mi, ITG *ncmat_, ITG *nstate_, double *cs, ITG *mcs, + ITG *nkon, double **enerp, double *xbounold, + double *xforcold, double *xloadold, + char *amname, double *amta, ITG *namta, + ITG *nam, ITG *iamforc, ITG *iamload, + ITG *iamt1, ITG *iamboun, double *ttime, char *output, + char *set, ITG *nset, ITG *istartset, + ITG *iendset, ITG *ialset, ITG *nprint, char *prlab, + char *prset, ITG *nener, double *trab, + ITG *inotr, ITG *ntrans, double *fmpc, ITG *ipobody, ITG *ibody, + double *xbody, ITG *nbody, double *xbodyold, double *timepar, + double *thicke, char *jobnamec, char *tieset, ITG *ntie, + ITG *istep, ITG *nmat, ITG *ielprop, double *prop, char *typeboun, + ITG *mortar, ITG *mpcinfo, double *tietol, ITG *ics, + char *orname, ITG *itempuser, double *t0g, double *t1g, + /* Adapter: Add variables for the participant name and the config file */ + char *preciceParticipantName, char *configFilename) +{ + + char description[13] = " ", *lakon = NULL, stiffmatrix[132] = "", + fneig[132] = "", jobnamef[396] = "", *labmpc2 = NULL; + + ITG *inum = NULL, k, *icol = NULL, *irow = NULL, ielas = 0, icmd = 0, iinc = 1, nasym = 0, i, j, ic, ir, + mass[2] = {0, 0}, stiffness = 1, buckling = 0, rhsi = 1, intscheme = 0, *ncocon = NULL, + *nshcon = NULL, mode = -1, noddiam = -1, coriolis = 0, iout, + *itg = NULL, ntg = 0, symmetryflag = 0, inputformat = 0, ngraph = 1, im, + mt = mi[1] + 1, ne0, *integerglob = NULL, iglob = 0, *ipneigh = NULL, *neigh = NULL, + icfd = 0, *inomat = NULL, *islavact = NULL, *islavnode = NULL, *nslavnode = NULL, + *islavsurf = NULL, nretain, *iretain = NULL, *noderetain = NULL, *ndirretain = NULL, + nmethodl, nintpoint, ifacecount, memmpc_, mpcfree, icascade, maxlenmpc, + ncont = 0, *itietri = NULL, *koncont = NULL, nslavs = 0, ismallsliding = 0, + *itiefac = NULL, *imastnode = NULL, *nmastnode = NULL, *imastop = NULL, iitsta, + *iponoels = NULL, *inoels = NULL, *ipe = NULL, *ime = NULL, iit = -1, iflagact = 0, + icutb = 0, *kon = NULL, *ipkon = NULL, *ielmat = NULL, ialeatoric = 0, kscale = 1, + *iponoel = NULL, *inoel = NULL, zero = 0, nherm = 1, nev = *nforc, node, idir, + *ielorien = NULL, network = 0, nrhs = 1, iperturbsav, mscalmethod = 0, *jqw = NULL, + *iroww = NULL, nzsw, *islavelinv = NULL, *irowtloc = NULL, *jqtloc = NULL, nboun2, + *ndirboun2 = NULL, *nodeboun2 = NULL, nmpc2, *ipompc2 = NULL, *nodempc2 = NULL, + *ikboun2 = NULL, *ilboun2 = NULL, *ikmpc2 = NULL, *ilmpc2 = NULL, mortartrafoflag = 0; + + double *stn = NULL, *v = NULL, *een = NULL, cam[5], *xstiff = NULL, *stiini = NULL, *tper, + *f = NULL, *fn = NULL, qa[4], *fext = NULL, *epn = NULL, *xstateini = NULL, + *vini = NULL, *stx = NULL, *enern = NULL, *xbounact = NULL, *xforcact = NULL, + *xloadact = NULL, *t1act = NULL, *ampli = NULL, *xstaten = NULL, *eei = NULL, + *enerini = NULL, *cocon = NULL, *shcon = NULL, *physcon = NULL, *qfx = NULL, + *qfn = NULL, sigma = 0., *cgr = NULL, *xbodyact = NULL, *vr = NULL, *vi = NULL, + *stnr = NULL, *stni = NULL, *vmax = NULL, *stnmax = NULL, *springarea = NULL, + *eenmax = NULL, *fnr = NULL, *fni = NULL, *emn = NULL, *clearini = NULL, ptime, + *emeini = NULL, *doubleglob = NULL, *au = NULL, *ad = NULL, *b = NULL, *aub = NULL, + *adb = NULL, *pslavsurf = NULL, *pmastsurf = NULL, *cdn = NULL, *cdnr = NULL, + *cdni = NULL, *submatrix = NULL, *xnoels = NULL, *cg = NULL, *straight = NULL, + *areaslav = NULL, *xmastnor = NULL, theta = 0., *ener = NULL, *xstate = NULL, + *fnext = NULL, *energyini = NULL, *energy = NULL, *d = NULL, alea = 0.1, *smscale = NULL, + *auw = NULL, *autloc = NULL, *xboun2 = NULL, *coefmpc2 = NULL; + + FILE *f1, *f2; + +#ifdef SGI + ITG token; +#endif + + /* dummy arguments for the results call */ + + double *veold = NULL, *accold = NULL, bet, gam, dtime, time, reltime = 1.; + + irow = *irowp; + ener = *enerp; + xstate = *xstatep; + ipkon = *ipkonp; + lakon = *lakonp; + kon = *konp; + ielmat = *ielmatp; + ielorien = *ielorienp; + icol = *icolp; + + for (k = 0; k < 3; k++) { + strcpy1(&jobnamef[k * 132], &jobnamec[k * 132], 132); + } + + tper = &timepar[1]; + + time = *tper; + dtime = *tper; + + ne0 = *ne; + + /* preCICE Adapter: Initialize the Calculix data structure */ + struct SimulationData simulationData = { + .ialset = ialset, + .ielmat = ielmat, + .istartset = istartset, + .iendset = iendset, + .kon = kon, + .ipkon = ipkon, + .lakon = lakon, + .co = co, + .set = set, + .nset = *nset, + .ikboun = ikboun, + // .ikforc = ikforc, + .ilboun = ilboun, + // .ilforc = ilforc, + .nboun = *nboun, + .nforc = *nforc, + .nelemload = nelemload, + .nload = *nload, + .sideload = sideload, + .mt = mt, + .nk = *nk, + .theta = &theta, + // .dtheta = &dtheta, + .tper = tper, + .nmethod = nmethod, + .xload = xload, + .xforc = xforc, + .xboun = xboun, + .ntmat_ = ntmat_, + .vold = vold, + .veold = veold, + .fn = fn, + .cocon = cocon, + .ncocon = ncocon, + .mi = mi, + .ne = *ne + }; + + /* preCICE Adapter: Initialize */ + printf("configFilename: %s\n", configFilename); + Precice_Setup(configFilename, preciceParticipantName, &simulationData); + Precice_AdjustSolverTimestep(&simulationData); + + /* determining the global values to be used as boundary conditions + for a submodel */ + + /* iglob=-1 if global results are from a *FREQUENCY calculation + iglob=0 if no global results are used by boundary conditions + iglob=1 if global results are from a *STATIC calculation */ + + ITG irefine = 0; + getglobalresults(&jobnamec[396], &integerglob, &doubleglob, nboun, iamboun, xboun, + nload, sideload, iamload, &iglob, nforc, iamforc, xforc, + ithermal, nk, t1, iamt1, &sigma, &irefine); + + /* reading temperatures from frd-file */ + + if ((itempuser[0] == 2) && (itempuser[1] != itempuser[2])) { + utempread(t1, &itempuser[2], jobnamec); + } + + /* allocating fields for the actual external loading */ + + NNEW(xbounact, double, *nboun); + for (k = 0; k < *nboun; ++k) { + xbounact[k] = xbounold[k]; + } + NNEW(xforcact, double, *nforc); + NNEW(xloadact, double, 2 * *nload); + NNEW(xbodyact, double, 7 * *nbody); + /* copying the rotation axis and/or acceleration vector */ + for (k = 0; k < 7 * *nbody; k++) { + xbodyact[k] = xbody[k]; + } + if (*ithermal == 1) { + NNEW(t1act, double, *nk); + for (k = 0; k < *nk; ++k) { + t1act[k] = t1old[k]; + } + } + + /* assigning the body forces to the elements */ + + /* if(*nbody>0){ + ifreebody=*ne+1; + NNEW(ipobody,ITG,2*ifreebody**nbody); + for(k=1;k<=*nbody;k++){ + FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset, + iendset,ialset,&inewton,nset,&ifreebody,&k)); + RENEW(ipobody,ITG,2*(*ne+ifreebody)); + } + RENEW(ipobody,ITG,2*(ifreebody-1)); + }*/ + + /* contact conditions */ + + // if(*icontact==1){ + if (*mortar > -2) { + + memmpc_ = mpcinfo[0]; + mpcfree = mpcinfo[1]; + icascade = mpcinfo[2]; + maxlenmpc = mpcinfo[3]; + + inicont(nk, &ncont, ntie, tieset, nset, set, istartset, iendset, ialset, &itietri, + lakon, ipkon, kon, &koncont, &nslavs, tietol, &ismallsliding, &itiefac, + &islavsurf, &islavnode, &imastnode, &nslavnode, &nmastnode, + mortar, &imastop, nkon, &iponoels, &inoels, &ipe, &ime, ne, &ifacecount, + iperturb, ikboun, nboun, co, istep, &xnoels); + + if (ncont != 0) { + + NNEW(cg, double, 3 * ncont); + NNEW(straight, double, 16 * ncont); + + /* 11 instead of 10: last position is reserved for the + local contact spring element number; needed as + pointer into springarea */ + + if (*mortar == 0) { + RENEW(kon, ITG, *nkon + 11 * nslavs); + NNEW(springarea, double, 2 * nslavs); + if (*nener == 1) { + RENEW(ener, double, 2 * mi[0] * (*ne + nslavs)); + DMEMSET(ener, 2 * mi[0] * *ne, 2 * mi[0] * (*ne + nslavs), 0.); + } + RENEW(ipkon, ITG, *ne + nslavs); + RENEW(lakon, char, 8 * (*ne + nslavs)); + + if (*norien > 0) { + RENEW(ielorien, ITG, mi[2] * (*ne + nslavs)); + for (k = mi[2] * *ne; k < mi[2] * (*ne + nslavs); k++) + ielorien[k] = 0; + } + + RENEW(ielmat, ITG, mi[2] * (*ne + nslavs)); + for (k = mi[2] * *ne; k < mi[2] * (*ne + nslavs); k++) + ielmat[k] = 1; + + if (nslavs != 0) { + RENEW(xstate, double, *nstate_ *mi[0] * (*ne + nslavs)); + for (k = *nstate_ * mi[0] * *ne; k < *nstate_ * mi[0] * (*ne + nslavs); k++) { + xstate[k] = 0.; + } + } + + NNEW(areaslav, double, ifacecount); + NNEW(xmastnor, double, 3 * nmastnode[*ntie]); + } else if (*mortar == 1) { + NNEW(islavact, ITG, nslavnode[*ntie]); + DMEMSET(islavact, 0, nslavnode[*ntie], 1); + NNEW(clearini, double, 3 * 9 * ifacecount); + NNEW(xmastnor, double, 3 * nmastnode[*ntie]); + + nintpoint = 0; + + precontact(&ncont, ntie, tieset, nset, set, istartset, + iendset, ialset, itietri, lakon, ipkon, kon, koncont, ne, + cg, straight, co, vold, istep, &iinc, &iit, itiefac, + islavsurf, islavnode, imastnode, nslavnode, nmastnode, + imastop, mi, ipe, ime, tietol, &iflagact, + &nintpoint, &pslavsurf, xmastnor, cs, mcs, ics, clearini, + &nslavs); + + /* changing the dimension of element-related fields */ + + RENEW(kon, ITG, *nkon + 22 * nintpoint); + RENEW(springarea, double, 2 * nintpoint); + RENEW(pmastsurf, double, 6 * nintpoint); + + if (*nener == 1) { + RENEW(ener, double, 2 * mi[0] * (*ne + nintpoint)); + DMEMSET(ener, 2 * mi[0] * *ne, 2 * mi[0] * (*ne + nintpoint), 0.); + } + RENEW(ipkon, ITG, *ne + nintpoint); + RENEW(lakon, char, 8 * (*ne + nintpoint)); + + if (*norien > 0) { + RENEW(ielorien, ITG, mi[2] * (*ne + nintpoint)); + for (k = mi[2] * *ne; k < mi[2] * (*ne + nintpoint); k++) + ielorien[k] = 0; + } + RENEW(ielmat, ITG, mi[2] * (*ne + nintpoint)); + for (k = mi[2] * *ne; k < mi[2] * (*ne + nintpoint); k++) + ielmat[k] = 1; + + /* interpolating the state variables */ + + if (*nstate_ != 0) { + + RENEW(xstate, double, *nstate_ *mi[0] * (ne0 + nintpoint)); + for (k = *nstate_ * mi[0] * ne0; k < *nstate_ * mi[0] * (ne0 + nintpoint); k++) { + xstate[k] = 0.; + } + + RENEW(xstateini, double, *nstate_ *mi[0] * (ne0 + nintpoint)); + for (k = 0; k < *nstate_ * mi[0] * (ne0 + nintpoint); ++k) { + xstateini[k] = xstate[k]; + } + } + } + + /* generating contact spring elements */ + + contact(&ncont, ntie, tieset, nset, set, istartset, iendset, + ialset, itietri, lakon, ipkon, kon, koncont, ne, cg, straight, nkon, + co, vold, ielmat, cs, elcon, istep, &iinc, &iit, ncmat_, ntmat_, + &ne0, nmethod, + iperturb, ikboun, nboun, mi, imastop, nslavnode, islavnode, islavsurf, + itiefac, areaslav, iponoels, inoels, springarea, tietol, &reltime, + imastnode, nmastnode, xmastnor, filab, mcs, ics, &nasym, + xnoels, mortar, pslavsurf, pmastsurf, clearini, &theta, + xstateini, xstate, nstate_, &icutb, &ialeatoric, jobnamef, + &alea, auw, jqw, iroww, &nzsw); + + printf("number of contact spring elements=%" ITGFORMAT "\n\n", *ne - ne0); + + /* determining the structure of the stiffness/mass matrix */ + + remastructar(ipompc, &coefmpc, &nodempc, nmpc, + &mpcfree, nodeboun, ndirboun, nboun, ikmpc, ilmpc, ikboun, ilboun, + labmpc, nk, &memmpc_, &icascade, &maxlenmpc, + kon, ipkon, lakon, ne, nactdof, icol, jq, &irow, isolver, + neq, nzs, nmethod, ithermal, iperturb, mass, mi, ics, cs, + mcs, mortar, typeboun, &iit, &network, iexpl); + } + + /* field for initial values of state variables (needed for contact */ + + if ((*nstate_ != 0) && ((*mortar == 0) || (ncont == 0))) { + NNEW(xstateini, double, *nstate_ *mi[0] * (ne0 + nslavs)); + for (k = 0; k < *nstate_ * mi[0] * (ne0 + nslavs); ++k) { + xstateini[k] = xstate[k]; + } + } + } + + /* allocating a field for the instantaneous amplitude */ + + NNEW(ampli, double, *nam); + + FORTRAN(tempload, (xforcold, xforc, xforcact, iamforc, nforc, xloadold, xload, + xloadact, iamload, nload, ibody, xbody, nbody, xbodyold, xbodyact, + t1old, t1, t1act, iamt1, nk, amta, + namta, nam, ampli, &time, &reltime, ttime, &dtime, ithermal, nmethod, + xbounold, xboun, xbounact, iamboun, nboun, + nodeboun, ndirboun, nodeforc, ndirforc, istep, &iinc, + co, vold, itg, &ntg, amname, ikboun, ilboun, nelemload, sideload, mi, + ntrans, trab, inotr, veold, integerglob, doubleglob, tieset, istartset, + iendset, ialset, ntie, nmpc, ipompc, ikmpc, ilmpc, nodempc, coefmpc, + ipobody, iponoel, inoel, ipkon, kon, ielprop, prop, ielmat, + shcon, nshcon, rhcon, nrhcon, cocon, ncocon, ntmat_, lakon, + set, nset)); + + /* determining the internal forces and the stiffness coefficients */ + + NNEW(f, double, *neq); + + /* allocating a field for the stiffness matrix */ + + NNEW(xstiff, double, (long long) 27 * mi[0] * *ne); + + /* for a *STATIC,PERTURBATION analysis with submodel boundary + conditions from a *FREQUENCY analysis iperturb[0]=1 has to be + temporarily set to iperturb[0]=0 in order for f to be calculated in + resultsini and subsequent results* routines */ + + if ((*nmethod == 1) && (iglob < 0) && (iperturb[0] > 0)) { + iperturbsav = iperturb[0]; + iperturb[0] = 0; + } + + iout = -1; + NNEW(v, double, mt **nk); + NNEW(fn, double, mt **nk); + NNEW(stx, double, 6 * mi[0] * *ne); + NNEW(inum, ITG, *nk); + NNEW(eei, double, 6 * mi[0] * *ne); + + results(co, nk, kon, ipkon, lakon, ne, v, stn, inum, stx, + elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, + ielorien, norien, orab, ntmat_, t0, t1act, ithermal, + prestr, iprestr, filab, eme, emn, een, iperturb, + f, fn, nactdof, &iout, qa, vold, b, nodeboun, + ndirboun, xbounact, nboun, ipompc, + nodempc, coefmpc, labmpc, nmpc, nmethod, cam, neq, veold, accold, + &bet, &gam, &dtime, &time, ttime, plicon, nplicon, plkcon, nplkcon, + xstateini, xstiff, xstate, npmat_, epn, matname, mi, &ielas, + &icmd, ncmat_, nstate_, stiini, vini, ikboun, ilboun, ener, enern, + emeini, xstaten, eei, enerini, cocon, ncocon, set, nset, istartset, + iendset, ialset, nprint, prlab, prset, qfx, qfn, trab, inotr, ntrans, + fmpc, nelemload, nload, ikmpc, ilmpc, istep, &iinc, springarea, + &reltime, &ne0, thicke, shcon, nshcon, + sideload, xloadact, xloadold, &icfd, inomat, pslavsurf, pmastsurf, + mortar, islavact, cdn, islavnode, nslavnode, ntie, clearini, + islavsurf, ielprop, prop, energyini, energy, &kscale, iponoel, + inoel, nener, orname, &network, ipobody, xbodyact, ibody, typeboun, + itiefac, tieset, smscale, &mscalmethod, nbody, t0g, t1g, + islavelinv, autloc, irowtloc, jqtloc, &nboun2, + ndirboun2, nodeboun2, xboun2, &nmpc2, ipompc2, nodempc2, coefmpc2, + labmpc2, ikboun2, ilboun2, ikmpc2, ilmpc2, &mortartrafoflag, + &intscheme); + + simulationData.xstiff = xstiff; + simulationData.eei = eei; + simulationData.stx = stx; + + if (Precice_IsCouplingOngoing()) { + printf("Write coupling data\n"); + Precice_WriteCouplingData(&simulationData); + + printf("Advancing the coupling\n"); + Precice_Advance(&simulationData); + + printf("Read the coupling\n"); + Precice_ReadCouplingData(&simulationData); + } + + simulationData.eei = NULL; + simulationData.stx = NULL; + simulationData.xstiff = NULL; + + SFREE(v); + SFREE(fn); + SFREE(stx); + SFREE(inum); + SFREE(eei); + iout = 1; + + if ((*nmethod == 1) && (iglob < 0) && (iperturb[0] > 0)) { + iperturb[0] = iperturbsav; + } + + /* determining the system matrix and the external forces */ + + NNEW(ad, double, *neq); + NNEW(fext, double, *neq); + + + if (*nmethod == 11) { + + /* determining the nodes and the degrees of freedom in those nodes + belonging to the substructure */ + + NNEW(iretain, ITG, *nk); + NNEW(noderetain, ITG, *nk); + NNEW(ndirretain, ITG, *nk); + nretain = 0; + + for (i = 0; i < *nboun; i++) { + if (strcmp1(&typeboun[i], "C") == 0) { + iretain[nretain] = i + 1; + noderetain[nretain] = nodeboun[i]; + ndirretain[nretain] = ndirboun[i]; + nretain++; + } + } + + /* nretain!=0: substructure application */ + + RENEW(iretain, ITG, nretain); + RENEW(noderetain, ITG, nretain); + RENEW(ndirretain, ITG, nretain); + + /* creating the right size au */ + + NNEW(au, double, nzs[2]); + rhsi = 0; + nmethodl = 2; + + } else { + + /* linear static calculation */ + + NNEW(au, double, *nzs); + nmethodl = *nmethod; + + /* if submodel calculation with a global model obtained by + a *FREQUENCY calculation: replace stiffness matrix K by + K-sigma*M */ + + if (iglob < 0) { + mass[0] = 1; + NNEW(adb, double, *neq); + NNEW(aub, double, nzs[1]); + } + } + + mafillsmmain(co, nk, kon, ipkon, lakon, ne, nodeboun, ndirboun, xbounact, nboun, + ipompc, nodempc, coefmpc, nmpc, nodeforc, ndirforc, xforcact, + nforc, nelemload, sideload, xloadact, nload, xbodyact, ipobody, + nbody, cgr, ad, au, fext, nactdof, icol, jq, irow, neq, nzl, &nmethodl, + ikmpc, ilmpc, ikboun, ilboun, + elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, + ielorien, norien, orab, ntmat_, + t0, t1act, ithermal, prestr, iprestr, vold, iperturb, sti, + nzs, stx, adb, aub, iexpl, plicon, nplicon, plkcon, nplkcon, + xstiff, npmat_, &dtime, matname, mi, + ncmat_, mass, &stiffness, &buckling, &rhsi, &intscheme, physcon, + shcon, nshcon, cocon, ncocon, ttime, &time, istep, &iinc, &coriolis, + ibody, xloadold, &reltime, veold, springarea, nstate_, + xstateini, xstate, thicke, integerglob, doubleglob, + tieset, istartset, iendset, ialset, ntie, &nasym, pslavsurf, + pmastsurf, mortar, clearini, ielprop, prop, &ne0, fnext, &kscale, + iponoel, inoel, &network, ntrans, inotr, trab, smscale, &mscalmethod, + set, nset, islavelinv, autloc, irowtloc, jqtloc, &mortartrafoflag); + + /* check for negative Jacobians */ + + if (nmethodl == 0) + *nmethod = 0; + + if (nasym == 1) { + RENEW(au, double, 2 * nzs[1]); + symmetryflag = 2; + inputformat = 1; + + mafillsmasmain(co, nk, kon, ipkon, lakon, ne, nodeboun, + ndirboun, xbounact, nboun, + ipompc, nodempc, coefmpc, nmpc, nodeforc, ndirforc, xforcact, + nforc, nelemload, sideload, xloadact, nload, xbodyact, ipobody, + nbody, cgr, ad, au, fext, nactdof, icol, jq, irow, neq, nzl, + nmethod, ikmpc, ilmpc, ikboun, ilboun, + elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, + ielmat, ielorien, norien, orab, ntmat_, + t0, t1act, ithermal, prestr, iprestr, vold, iperturb, sti, + nzs, stx, adb, aub, iexpl, plicon, nplicon, plkcon, nplkcon, + xstiff, npmat_, &dtime, matname, mi, + ncmat_, mass, &stiffness, &buckling, &rhsi, &intscheme, + physcon, shcon, nshcon, cocon, ncocon, ttime, &time, istep, &iinc, + &coriolis, ibody, xloadold, &reltime, veold, springarea, nstate_, + xstateini, xstate, thicke, + integerglob, doubleglob, tieset, istartset, iendset, + ialset, ntie, &nasym, pslavsurf, pmastsurf, mortar, clearini, + ielprop, prop, &ne0, &kscale, iponoel, inoel, &network, set, nset); + } + + /* determining the right hand side */ + + NNEW(b, double, *neq); + for (k = 0; k < *neq; ++k) { + b[k] = fext[k] - f[k]; + } + SFREE(fext); + SFREE(f); + + /* generation of a substructure stiffness matrix */ + + if (*nmethod == 11) { + + /* factorizing the matrix */ + + if (*neq > 0) { + if (*isolver == 0) { +#ifdef SPOOLES + spooles_factor(ad, au, adb, aub, &sigma, icol, irow, neq, nzs, &symmetryflag, + &inputformat, &nzs[2]); +#else + printf(" *ERROR in linstatic: the SPOOLES library is not linked\n\n"); + FORTRAN(stop, ()); +#endif + } else if (*isolver == 7) { +#ifdef PARDISO + pardiso_factor(ad, au, adb, aub, &sigma, icol, irow, neq, nzs, + &symmetryflag, &inputformat, jq, &nzs[2]); +#else + printf(" *ERROR in linstatic: the PARDISO library is not linked\n\n"); + FORTRAN(stop, ()); +#endif + } else if (*isolver == 8) { +#ifdef PASTIX + pastix_factor_main(ad, au, adb, aub, &sigma, icol, irow, neq, nzs, + &symmetryflag, &inputformat, jq, &nzs[2]); +#else + printf(" *ERROR in linstatic: the PASTIX library is not linked\n\n"); + FORTRAN(stop, ()); +#endif + } + } + + /* solving the system of equations with appropriate rhs */ + + /* substructure calculations */ + + NNEW(submatrix, double, nretain *nretain); + + for (i = 0; i < nretain; i++) { + DMEMSET(b, 0, *neq, 0.); + ic = *neq + iretain[i] - 1; + for (j = jq[ic] - 1; j < jq[ic + 1] - 1; j++) { + ir = irow[j] - 1; + b[ir] -= au[j]; + } + + /* solving the system */ + + if (*neq > 0) { + if (*isolver == 0) { +#ifdef SPOOLES + spooles_solve(b, neq); +#endif + } else if (*isolver == 7) { +#ifdef PARDISO + pardiso_solve(b, neq, &symmetryflag, &inputformat, &nrhs); +#endif + + } else if (*isolver == 8) { +#ifdef PASTIX + pastix_solve(b, neq, &symmetryflag, &nrhs); +#endif + } + } + + /* calculating the internal forces */ + + NNEW(v, double, mt **nk); + NNEW(fn, double, mt **nk); + NNEW(stn, double, 6 * *nk); + NNEW(inum, ITG, *nk); + NNEW(stx, double, 6 * mi[0] * *ne); + + if (strcmp1(&filab[261], "E ") == 0) + NNEW(een, double, 6 * *nk); + if (strcmp1(&filab[2697], "ME ") == 0) + NNEW(emn, double, 6 * *nk); + if (strcmp1(&filab[522], "ENER") == 0) + NNEW(enern, double, *nk); + + NNEW(eei, double, 6 * mi[0] * *ne); + if (*nener == 1) { + NNEW(stiini, double, 6 * mi[0] * *ne); + NNEW(emeini, double, 6 * mi[0] * *ne); + NNEW(enerini, double, 2 * mi[0] * *ne); + } + + /* replacing the appropriate boundary value by unity */ + + xbounact[iretain[i] - 1] = 1.; + + results(co, nk, kon, ipkon, lakon, ne, v, stn, inum, stx, + elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, + ielorien, norien, orab, ntmat_, t0, t1act, ithermal, + prestr, iprestr, filab, eme, emn, een, iperturb, + f, fn, nactdof, &iout, qa, vold, b, nodeboun, ndirboun, + xbounact, nboun, ipompc, + nodempc, coefmpc, labmpc, nmpc, nmethod, cam, neq, veold, + accold, &bet, + &gam, &dtime, &time, ttime, plicon, nplicon, plkcon, nplkcon, + xstateini, xstiff, xstate, npmat_, epn, matname, mi, &ielas, &icmd, + ncmat_, nstate_, stiini, vini, ikboun, ilboun, ener, enern, emeini, + xstaten, eei, enerini, cocon, ncocon, set, nset, istartset, iendset, + ialset, nprint, prlab, prset, qfx, qfn, trab, inotr, ntrans, fmpc, + nelemload, nload, ikmpc, ilmpc, istep, &iinc, springarea, &reltime, + &ne0, thicke, shcon, nshcon, + sideload, xloadact, xloadold, &icfd, inomat, pslavsurf, pmastsurf, + mortar, islavact, cdn, islavnode, nslavnode, ntie, clearini, + islavsurf, ielprop, prop, energyini, energy, &kscale, iponoel, + inoel, nener, orname, &network, ipobody, xbodyact, ibody, typeboun, + itiefac, tieset, smscale, &mscalmethod, nbody, t0g, t1g, + islavelinv, autloc, irowtloc, jqtloc, &nboun2, + ndirboun2, nodeboun2, xboun2, &nmpc2, ipompc2, nodempc2, coefmpc2, + labmpc2, ikboun2, ilboun2, ikmpc2, ilmpc2, &mortartrafoflag, + &intscheme); + + xbounact[iretain[i] - 1] = 0.; + + SFREE(v); + SFREE(stn); + SFREE(inum); + SFREE(stx); + SFREE(eei); + + if (strcmp1(&filab[261], "E ") == 0) + SFREE(een); + if (strcmp1(&filab[2697], "ME ") == 0) + SFREE(emn); + if (strcmp1(&filab[522], "ENER") == 0) + SFREE(enern); + + SFREE(eei); + if (*nener == 1) { + SFREE(stiini); + SFREE(emeini); + SFREE(enerini); + } + + /* storing the internal forces in the substructure + stiffness matrix */ + + for (j = 0; j < nretain; j++) { + submatrix[i * nretain + j] = fn[mt * (noderetain[j] - 1) + ndirretain[j]]; + } + + SFREE(fn); + } + + /* preCICE cleanup*/ + precicec_finalize(); + + /* cleanup */ + + if (*neq > 0) { + if (*isolver == 0) { +#ifdef SPOOLES + spooles_cleanup(); +#endif + } else if (*isolver == 7) { +#ifdef PARDISO + pardiso_cleanup(&neq[0], &symmetryflag, &inputformat); +#endif + } else if (*isolver == 8) { +#ifdef PASTIX +#endif + } + } + + SFREE(iretain); + + FORTRAN(writesubmatrix, (submatrix, noderetain, ndirretain, &nretain, jobnamec)); + + SFREE(submatrix); + SFREE(noderetain); + SFREE(ndirretain); + + SFREE(au); + SFREE(ad); + SFREE(b); + + SFREE(xbounact); + SFREE(xforcact); + SFREE(xloadact); + SFREE(t1act); + SFREE(ampli); + SFREE(xbodyact); + + // if(*nbody>0) SFREE(ipobody); + + SFREE(xstiff); + + if (iglob != 0) { + SFREE(integerglob); + SFREE(doubleglob); + } + + return; + + } else if (*nmethod != 0) { + /* linear static applications */ + printf("Solve the linear static problem\n"); + + if (*isolver == 0) { +#ifdef SPOOLES + spooles(ad, au, adb, aub, &sigma, b, icol, irow, neq, nzs, &symmetryflag, + &inputformat, &nzs[2]); + printf("After spooles\n"); +#else + printf(" *ERROR in linstatic: the SPOOLES library is not linked\n\n"); + FORTRAN(stop, ()); +#endif + } else if ((*isolver == 2) || (*isolver == 3)) { + if (nasym > 0) { + printf(" *ERROR in nonlingeo: the iterative solver cannot be used for asymmetric matrices\n\n"); + FORTRAN(stop, ()); + } + preiter(ad, &au, b, &icol, &irow, neq, nzs, isolver, iperturb); + } else if (*isolver == 4) { +#ifdef SGI + if (nasym > 0) { + printf(" *ERROR in nonlingeo: the SGI solver cannot be used for asymmetric matrices\n\n"); + FORTRAN(stop, ()); + } + token = 1; + sgi_main(ad, au, adb, aub, &sigma, b, icol, irow, neq, nzs, token); +#else + printf(" *ERROR in linstatic: the SGI library is not linked\n\n"); + FORTRAN(stop, ()); +#endif + } else if (*isolver == 5) { +#ifdef TAUCS + if (nasym > 0) { + printf(" *ERROR in nonlingeo: the TAUCS solver cannot be used for asymmetric matrices\n\n"); + FORTRAN(stop, ()); + } + tau(ad, &au, adb, aub, &sigma, b, icol, &irow, neq, nzs); +#else + printf(" *ERROR in linstatic: the TAUCS library is not linked\n\n"); + FORTRAN(stop, ()); +#endif + } else if (*isolver == 7) { +#ifdef PARDISO + pardiso_main(ad, au, adb, aub, &sigma, b, icol, irow, neq, nzs, + &symmetryflag, &inputformat, jq, &nzs[2], &nrhs); +#else + printf(" *ERROR in linstatic: the PARDISO library is not linked\n\n"); + FORTRAN(stop, ()); +#endif + } else if (*isolver == 8) { +#ifdef PASTIX + pastix_main(ad, au, adb, aub, &sigma, b, icol, irow, neq, nzs, + &symmetryflag, &inputformat, jq, &nzs[2], &nrhs); +#else + printf(" *ERROR in linstatic: the PASTIX library is not linked\n\n"); + FORTRAN(stop, ()); +#endif + } + + /* saving of ad and au for sensitivity analysis */ + + for (i = 0; i < *ntie; i++) { + if (strcmp1(&tieset[i * 243 + 80], "D") == 0) { + + // strcpy2(stiffmatrix, jobnamec, 132); + strcat(stiffmatrix, ".stm"); + + if ((f1 = fopen(stiffmatrix, "wb")) == NULL) { + printf(" *ERROR in linstatic: cannot open stiffness matrix file for writing..."); + exit(0); + } + + /* storing the stiffness matrix */ + printf("DEBUGGING STEP 2\n"); + + /* nzs,irow,jq and icol have to be stored too, since the static analysis + can involve contact, whereas in the sensitivity analysis contact is not + taken into account while determining the structure of the stiffness + matrix (in mastruct.c) + */ + + if (fwrite(&nasym, sizeof(ITG), 1, f1) != 1) { + printf(" *ERROR saving the symmetry flag to the stiffness matrix file..."); + exit(0); + } + if (fwrite(nzs, sizeof(ITG), 3, f1) != 3) { + printf(" *ERROR saving the number of subdiagonal nonzeros to the stiffness matrix file..."); + exit(0); + } + if (fwrite(irow, sizeof(ITG), nzs[2], f1) != nzs[2]) { + printf(" *ERROR saving irow to the stiffness matrix file..."); + exit(0); + } + if (fwrite(jq, sizeof(ITG), neq[1] + 1, f1) != neq[1] + 1) { + printf(" *ERROR saving jq to the stiffness matrix file..."); + exit(0); + } + if (fwrite(icol, sizeof(ITG), neq[1], f1) != neq[1]) { + printf(" *ERROR saving icol to the stiffness matrix file..."); + exit(0); + } + if (fwrite(ad, sizeof(double), neq[1], f1) != neq[1]) { + printf(" *ERROR saving the diagonal of the stiffness matrix to the stiffness matrix file..."); + exit(0); + } + if (fwrite(au, sizeof(double), nzs[2], f1) != nzs[2]) { + printf(" *ERROR saving the off-diagonal terms of the stiffness matrix to the tiffness matrix file..."); + exit(0); + } + fclose(f1); + + break; + } + } + + SFREE(ad); + SFREE(au); + if (iglob < 0) { + SFREE(adb); + SFREE(aub); + } + + /* calculating the displacements and the stresses and storing */ + /* the results in frd format for each valid eigenmode */ + + NNEW(v, double, mt **nk); + NNEW(fn, double, mt **nk); + NNEW(stn, double, 6 * *nk); + NNEW(inum, ITG, *nk); + NNEW(stx, double, 6 * mi[0] * *ne); + + if (strcmp1(&filab[261], "E ") == 0) + NNEW(een, double, 6 * *nk); + if (strcmp1(&filab[2697], "ME ") == 0) + NNEW(emn, double, 6 * *nk); + if (strcmp1(&filab[522], "ENER") == 0) + NNEW(enern, double, *nk); + if (strcmp1(&filab[2175], "CONT") == 0) + NNEW(cdn, double, 6 * *nk); + + NNEW(eei, double, 6 * mi[0] * *ne); + if (*nener == 1) { + NNEW(stiini, double, 6 * mi[0] * *ne); + NNEW(emeini, double, 6 * mi[0] * *ne); + NNEW(enerini, double, 2 * mi[0] * *ne); + } + + results(co, nk, kon, ipkon, lakon, ne, v, stn, inum, stx, + elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, + ielorien, norien, orab, ntmat_, t0, t1act, ithermal, + prestr, iprestr, filab, eme, emn, een, iperturb, + f, fn, nactdof, &iout, qa, vold, b, nodeboun, ndirboun, xbounact, nboun, + ipompc, + nodempc, coefmpc, labmpc, nmpc, nmethod, cam, neq, veold, accold, &bet, + &gam, &dtime, &time, ttime, plicon, nplicon, plkcon, nplkcon, + xstateini, xstiff, xstate, npmat_, epn, matname, mi, &ielas, &icmd, + ncmat_, nstate_, stiini, vini, ikboun, ilboun, ener, enern, emeini, + xstaten, eei, enerini, cocon, ncocon, set, nset, istartset, iendset, + ialset, nprint, prlab, prset, qfx, qfn, trab, inotr, ntrans, fmpc, + nelemload, nload, ikmpc, ilmpc, istep, &iinc, springarea, &reltime, + &ne0, thicke, shcon, nshcon, + sideload, xloadact, xloadold, &icfd, inomat, pslavsurf, pmastsurf, + mortar, islavact, cdn, islavnode, nslavnode, ntie, clearini, + islavsurf, ielprop, prop, energyini, energy, &kscale, iponoel, + inoel, nener, orname, &network, ipobody, xbodyact, ibody, typeboun, + itiefac, tieset, smscale, &mscalmethod, nbody, t0g, t1g, + islavelinv, autloc, irowtloc, jqtloc, &nboun2, + ndirboun2, nodeboun2, xboun2, &nmpc2, ipompc2, nodempc2, coefmpc2, + labmpc2, ikboun2, ilboun2, ikmpc2, ilmpc2, &mortartrafoflag, + &intscheme); + + SFREE(eei); + if (*nener == 1) { + SFREE(stiini); + SFREE(emeini); + SFREE(enerini); + } + + simulationData.xstiff = xstiff; + simulationData.eei = eei; + simulationData.stx = stx; + + if (Precice_IsCouplingOngoing()) { + printf("Write coupling data\n"); + Precice_WriteCouplingData(&simulationData); + + printf("Advance coupling\n"); + Precice_Advance(&simulationData); + + // Last read command is not necessary as the data read is never used. + //printf("Read coupling data\n"); + //Precice_ReadCouplingData(&simulationData); + } + + simulationData.eei = NULL; + simulationData.stx = NULL; + simulationData.xstiff = NULL; + + memcpy(&vold[0], &v[0], sizeof(double) * mt * *nk); + memcpy(&sti[0], &stx[0], sizeof(double) * 6 * mi[0] * ne0); + + ++*kode; + + /* for cyclic symmetric sectors: duplicating the results */ + + if (*mcs > 0) { + ptime = *ttime + time; + frdcyc(co, nk, kon, ipkon, lakon, ne, v, stn, inum, nmethod, kode, filab, een, t1act, + fn, &ptime, epn, ielmat, matname, cs, mcs, nkon, enern, xstaten, + nstate_, istep, &iinc, iperturb, ener, mi, output, ithermal, + qfn, ialset, istartset, iendset, trab, inotr, ntrans, orab, + ielorien, norien, sti, veold, &noddiam, set, nset, emn, thicke, + jobnamec, &ne0, cdn, mortar, nmat, qfx, ielprop, prop); + } else { + if (strcmp1(&filab[1044], "ZZS") == 0) { + NNEW(neigh, ITG, 40 * *ne); + NNEW(ipneigh, ITG, *nk); + } + ptime = *ttime + time; + frd(co, nk, kon, ipkon, lakon, ne, v, stn, inum, nmethod, + kode, filab, een, t1act, fn, &ptime, epn, ielmat, matname, enern, xstaten, + nstate_, istep, &iinc, ithermal, qfn, &mode, &noddiam, trab, inotr, + ntrans, orab, ielorien, norien, description, ipneigh, neigh, + mi, stx, vr, vi, stnr, stni, vmax, stnmax, &ngraph, veold, ener, ne, + cs, set, nset, istartset, iendset, ialset, eenmax, fnr, fni, emn, + thicke, jobnamec, output, qfx, cdn, mortar, cdnr, cdni, nmat, ielprop, + prop, sti); + if (strcmp1(&filab[1044], "ZZS") == 0) { + SFREE(ipneigh); + SFREE(neigh); + } + } + + /* updating the .sta file */ + + iitsta = 1; + FORTRAN(writesta, (istep, &iinc, &icutb, &iitsta, ttime, &time, &dtime)); + + SFREE(v); + SFREE(stn); + SFREE(inum); + SFREE(b); + SFREE(stx); + SFREE(fn); + + if (strcmp1(&filab[261], "E ") == 0) + SFREE(een); + if (strcmp1(&filab[2697], "ME ") == 0) + SFREE(emn); + if (strcmp1(&filab[522], "ENER") == 0) + SFREE(enern); + if (strcmp1(&filab[2175], "CONT") == 0) + SFREE(cdn); + + } else { + + /* error occurred in mafill: storing the geometry in frd format */ + + ++*kode; + NNEW(inum, ITG, *nk); + for (k = 0; k < *nk; k++) + inum[k] = 1; + if (strcmp1(&filab[1044], "ZZS") == 0) { + NNEW(neigh, ITG, 40 * *ne); + NNEW(ipneigh, ITG, *nk); + } + ptime = *ttime + time; + frd(co, nk, kon, ipkon, lakon, ne, v, stn, inum, nmethod, + kode, filab, een, t1, fn, &ptime, epn, ielmat, matname, enern, xstaten, + nstate_, istep, &iinc, ithermal, qfn, &mode, &noddiam, trab, inotr, + ntrans, orab, ielorien, norien, description, ipneigh, neigh, + mi, sti, vr, vi, stnr, stni, vmax, stnmax, &ngraph, veold, ener, ne, + cs, set, nset, istartset, iendset, ialset, eenmax, fnr, fni, emn, + thicke, jobnamec, output, qfx, cdn, mortar, cdnr, cdni, nmat, ielprop, + prop, sti); + if (strcmp1(&filab[1044], "ZZS") == 0) { + SFREE(ipneigh); + SFREE(neigh); + } + SFREE(inum); + FORTRAN(stop, ()); + } + + if (*mortar > -2) { + if (ncont != 0) { + *ne = ne0; + if (*nener == 1) + RENEW(ener, double, 2 * mi[0] * *ne); + RENEW(ipkon, ITG, *ne); + RENEW(lakon, char, 8 * *ne); + RENEW(kon, ITG, *nkon); + if (*norien > 0) { + RENEW(ielorien, ITG, mi[2] * *ne); + } + RENEW(ielmat, ITG, mi[2] * *ne); + SFREE(cg); + SFREE(straight); + SFREE(imastop); + SFREE(itiefac); + SFREE(islavnode); + SFREE(islavsurf); + SFREE(nslavnode); + SFREE(iponoels); + SFREE(inoels); + SFREE(imastnode); + SFREE(nmastnode); + SFREE(itietri); + SFREE(koncont); + SFREE(xnoels); + SFREE(springarea); + SFREE(xmastnor); + + if (*mortar == 0) { + SFREE(areaslav); + } else if (*mortar == 1) { + SFREE(pmastsurf); + SFREE(ipe); + SFREE(ime); + SFREE(pslavsurf); + SFREE(islavact); + SFREE(clearini); + } + } + mpcinfo[0] = memmpc_; + mpcinfo[1] = mpcfree; + mpcinfo[2] = icascade; + mpcinfo[3] = maxlenmpc; + } + + /* updating the loading at the end of the step; + important in case the amplitude at the end of the step + is not equal to one */ + + for (k = 0; k < *nboun; ++k) { + xbounold[k] = xbounact[k]; + } + for (k = 0; k < *nforc; ++k) { + xforcold[k] = xforcact[k]; + } + for (k = 0; k < 2 * *nload; ++k) { + xloadold[k] = xloadact[k]; + } + for (k = 0; k < 7 * *nbody; k = k + 7) { + xbodyold[k] = xbodyact[k]; + } + if (*ithermal == 1) { + for (k = 0; k < *nk; ++k) { + t1old[k] = t1act[k]; + } + for (k = 0; k < *nk; ++k) { + vold[mt * k] = t1act[k]; + } + } + + SFREE(xbounact); + SFREE(xforcact); + SFREE(xloadact); + SFREE(t1act); + SFREE(ampli); + SFREE(xbodyact); + + // if(*nbody>0) SFREE(ipobody); + + if (iglob != 0) { + SFREE(integerglob); + SFREE(doubleglob); + } + + *irowp = irow; + *enerp = ener; + *xstatep = xstate; + *ipkonp = ipkon; + *lakonp = lakon; + *konp = kon; + *ielmatp = ielmat; + *ielorienp = ielorien; + *icolp = icol; + + (*ttime) += (*tper); + + /* preCICE Adapter: Free the memory */ + Precice_FreeData(&simulationData); + + return; +}