Skip to content

Commit

Permalink
Trying to get gauss point coorindates from the FORTRAN side
Browse files Browse the repository at this point in the history
  • Loading branch information
@IshaanDesai
IshaanDesai committed Aug 19, 2024
1 parent a56a7cf commit 38fa1b8
Show file tree
Hide file tree
Showing 4 changed files with 116 additions and 100 deletions.
41 changes: 9 additions & 32 deletions adapter/CCXHelpers.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,11 +93,15 @@ void getNodeCoordinates(ITG *nodes, ITG numNodes, int dim, double *co, double *v
}
}

//, d
void getElementGaussPointCoordinates(int numElements, int numGPTotal, int *elementIDs, double *co, ITG *kon, char *lakon, ITG *ipkon, int *gp_id, double *gp_coord)
{
// FORTRAN(getelementgausspointcoords, (&numElements, &numGPTotal, elementIDs, co, &lakon, kon, ipkon, gp_id, gp_coord));
}
// void getElementGaussPointCoordinates(int numElements, int numGPTotal, int *elementIDs, double *co, ITG *kon, char *lakon, ITG *ipkon, int *gp_id, double *gp_coord)
// {
// FORTRAN(getelementgausspointcoords, (&numElements, &numGPTotal, elementIDs, co, &lakon, kon, ipkon, gp_id, gp_coord));

// printf("Element coordinates inside CCXHelpers.c\n");
// for (int j = 0; j < numGPTotal; j++) {
// printf(" %d, element coordinates: %f, %f, %f \n", j, gp_coord[j * 3], gp_coord[j * 3 + 1], gp_coord[j * 3 + 2]);
// }
// }

void getNodeTemperatures(ITG *nodes, ITG numNodes, double *v, int mt, double *temperatures)
{
Expand Down Expand Up @@ -502,33 +506,6 @@ void setNodeDisplacements(double *displacements, ITG numNodes, int dim, int *xbo
}
}

void setElementXstiff(int nelem, ITG *mi, double *cmatData, double *xstiff)
{
printf("bbbefore setting xstiff\n");

// int i, count, j,xstiffSize, nSize;
// xstiffSize = 27;
// nSize = mi[0]*nelem;
// cidx = 15;

printf("before setting xstiff\n");
for (int i = 0; i < mi[0] * nelem * 27; i++) {
xstiff[i] = 1.0;
}
printf("after setting xstiff\n");

// // Loop through all element and respective gauss points
// count=0;
// for (i = 0; i < nSize; i++) {
// j = i*xstiffSize+cidx;
// printf("idx: %ld\n",j);
// xstiff[j] = 1.0; //cmatData[count];
// xstiff[j+1] = 1.0; //cmatData[count+1];
// xstiff[j+2] = 1.0; //cmatData[count+2];
// count = count + 3;
// }
}

bool isSteadyStateSimulation(ITG *nmethod)
{
return *nmethod == 1;
Expand Down
56 changes: 40 additions & 16 deletions adapter/PreciceInterface.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -651,11 +651,11 @@ void PreciceInterface_Create(PreciceInterface *interface, SimulationData *sim, I
// Initialize pointers as NULL
interface->elementIDs = NULL;
interface->faceIDs = NULL;
interface->nodeIDs = NULL;
interface->faceCenterCoordinates = NULL;
interface->preciceFaceCenterIDs = NULL;
interface->nodeCoordinates = NULL;
interface->node2DCoordinates = NULL;
interface->nodeIDs = NULL;
interface->mapping2D3D = NULL;
interface->preciceNodeIDs = NULL;
interface->nodeScalarData = NULL;
Expand Down Expand Up @@ -767,25 +767,40 @@ void PreciceInterface_ConfigureElementsMesh(PreciceInterface *interface, Simulat
interface->elementIDs = malloc(interface->numElements * sizeof(ITG));
getElementsIDs(interface->elementSetID, sim->ialset, sim->istartset, sim->iendset, interface->elementIDs);

for (int j = 0; j < interface->numElements; j++) {
printf(" %d, element id: %d \n", j, interface->elementIDs[j]);
}

// Find guass point coordinates of the element -> Serves as mesh for data transfer
int numElt = interface->numElements;
interface->numIPTotal = 8 * interface->numElements; // Gauss point mesh coordinate -Each element 8 gauss points
interface->numIPTotal = 8 * interface->numElements; // Gauss point mesh coordinate - each element 8 gauss points
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;
interface->elemIPCoordinates[j * 3] = j;
interface->elemIPCoordinates[j * 3 + 1] = 0.0;
interface->elemIPCoordinates[j * 3 + 2] = 0.0;
interface->elemIPID[j] = j;
}

// getElementGaussPointCoordinates(interface->numElements, interface->numIPTotal, interface->elementIDs, sim->co,
// sim->kon, sim->lakon, sim->ipkon, interface->elemIPID, interface->elemIPCoordinates);
//getElementGaussPointCoordinates(interface->numElements, interface->numIPTotal, interface->elementIDs, sim->co,
// sim->kon, sim->lakon, sim->ipkon, interface->elemIPID, interface->elemIPCoordinates);

int *elem_ids = malloc(interface->numElements * sizeof(ITG));

printf("num elements: %d\n", interface->numElements);

int numElements = interface->numElements;

FORTRAN(getelementgausspointcoords, (&numElements,
interface->numIPTotal,
interface->elementIDs,
sim->co,
sim->lakon,
sim->kon,
sim->ipkon,
elem_ids,
interface->elemIPCoordinates));

free(elem_ids);

// print element coordinates
for (int j = 0; j < interface->numIPTotal; j++) {
printf(" %d, element coordinates: %f, %f, %f \n", j, interface->elemIPCoordinates[j * 3], interface->elemIPCoordinates[j * 3 + 1], interface->elemIPCoordinates[j * 3 + 2]);
}

precicec_setMeshVertices(interface->elementMeshName, interface->numIPTotal, interface->elemIPCoordinates, interface->elemIPID);
}
Expand Down Expand Up @@ -1127,9 +1142,11 @@ 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);
free(preciceInterface->node2DCoordinates);
free(preciceInterface->preciceNodeIDs);
free(preciceInterface->nodeScalarData);
free(preciceInterface->node2DScalarData);
Expand All @@ -1146,12 +1163,15 @@ void PreciceInterface_FreeData(PreciceInterface *preciceInterface)
free(preciceInterface->elementIPScalarData);
free(preciceInterface->elementIPVectorData);

// Quasi 2D-3D coupling
free(preciceInterface->mapping2D3D);
freeMapping(preciceInterface->mappingQuasi2D3D);

// Mesh names
free(preciceInterface->faceCentersMeshName);
free(preciceInterface->nodesMeshName);
free(preciceInterface->elementMeshName);
free(preciceInterface->couplingMeshName);

// Data names
free(preciceInterface->displacementDeltas);
Expand All @@ -1177,6 +1197,10 @@ void PreciceInterface_FreeData(PreciceInterface *preciceInterface)
free(preciceInterface->materialTangent5Data);
free(preciceInterface->materialTangent6Data);
free(preciceInterface->materialTangent7Data);
free(preciceInterface->rveid);
free(preciceInterface->modid);
free(preciceInterface->rucsize);
free(preciceInterface->conv);
}

void PreciceInterface_MultiscaleCheckpoint(SimulationData *sim)
Expand All @@ -1187,10 +1211,10 @@ void PreciceInterface_MultiscaleCheckpoint(SimulationData *sim)
// Write strain data
Precice_WriteCouplingData(sim);

// Advance time
Precice_Advance(sim);

// Read stress, material tangent data
Precice_ReadCouplingData(sim);

// Advance time
Precice_Advance(sim);
}
}
62 changes: 30 additions & 32 deletions getelementgausspointcoords.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -242,35 +242,39 @@ subroutine getelementgausspointcoords(NELEM, NGP, ELEM_IDS,CO,

implicit none


!Input/Output variables
integer :: NELEM ! Number of elements
integer :: NGP ! Total number of elements (nelem * 8)
integer :: ELEM_IDS(*) ! Element Ids - Need to increment by 1 for C->Fortran conversion
real(8) :: CO(3,*) ! Nodal coordinates of all nodes
integer :: NELEM ! Number of elements
integer :: NGP ! Total number of elements (nelem * 8)
integer :: ELEM_IDS(*) ! Element IDs
real(8) :: CO(3,*) ! Nodal coordinates of all nodes
character(8):: LAKON(*)
integer :: KON(*)
integer :: IPKON(*)

integer :: ELEM_GP_ID(*) ! GP ID
real :: ELEM_GP_COORD(*) ! GP coords

integer :: ELEM_GP_ID(*) ! GP ID
real :: ELEM_GP_COORD(*) ! GP coords

! !Internal variables
INTEGER(4) :: IEL, EL_ID, INDEXE, I,J,K,L, KONL(8), GP_ID
INTEGER(4) :: GP_TOT,IFLAG
CHARACTER(8) :: LAKONL
REAL(8) :: XL(3,8), XI, ET, ZE, xsj,shp(4,20)


include "gauss.f"

data iflag /1/

WRITE(*,*) "STEP 1"

WRITE(*,*) "NELEM: ", NELEM

GP_TOT = 1
DO IEL = 1, NELEM
WRITE(*,*) "IEL: ", IEL
EL_ID = ELEM_IDS(IEL)

WRITE(*,*) "EL_ID: ", EL_ID

LAKONL=LAKON(EL_ID)
INDEXE=ipkon(EL_ID)

Expand All @@ -286,7 +290,7 @@ subroutine getelementgausspointcoords(NELEM, NGP, ELEM_IDS,CO,
enddo
enddo

! Loop through gauss points fo each element
! Loop through gauss points of each element
DO GP_ID = 1,8
ELEM_GP_ID(GP_TOT) = GP_TOT-1

Expand All @@ -300,37 +304,31 @@ subroutine getelementgausspointcoords(NELEM, NGP, ELEM_IDS,CO,

! GAUSS POINT COORDINATES

! do k=1,3
! do l=1,8
! ELEM_GP_COORD(k,GP_TOT)=ELEM_GP_COORD(k,GP_TOT)
! & +xl(k,l)*shp(4,l)
! enddo
! enddo
!do k=1,3
! do l=1,8
! ELEM_GP_COORD(k*l)=ELEM_GP_COORD(k*l)
! & +xl(k,l)*shp(4,l)
! enddo
!enddo
ELEM_GP_COORD((GP_TOT-1)*3+1) = GP_TOT-1
! ELEM_GP_COORD(1,GP_TOT)=GP_TOT-1
write(*,*) ELEM_GP_ID(GP_TOT), ELEM_GP_COORD((GP_TOT-1)*3+1)

GP_TOT = GP_TOT+1



ENDDO










ENDDO

write(*,*) "ELEMENT COORDS IN FORTRAN"
write(*,*) "========================="
GP_TOT = 1
do iel = 1, NELEM
do j = 1, 8
write(*,*) ELEM_GP_ID(GP_TOT), ELEM_GP_COORD((GP_TOT-1)*3+1)
GP_TOT = GP_TOT+1
enddo
enddo

RETURN

end subroutine getelementgausspointcoords




Loading

0 comments on commit 38fa1b8

Please sign in to comment.