Skip to content

Commit

Permalink
added strain energy calculations
Browse files Browse the repository at this point in the history
  • Loading branch information
@chennachaos
chennachaos committed Nov 29, 2023
1 parent 9ef688d commit c52d701
Show file tree
Hide file tree
Showing 8 changed files with 186 additions and 65 deletions.
1 change: 1 addition & 0 deletions project/sampleinputs/IbeamVertical
View file
Original file line number Diff line number Diff line change
Expand Up @@ -859,6 +859,7 @@ immersed body output
3 336 2
4 336 1
4 336 2
100 1 1


control parameters
Expand Down
5 changes: 3 additions & 2 deletions src/HBsplines/HBSplineBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,13 +66,14 @@ class HBSplineBase: public Domain

vector<int> VacantBFs, nodes2divide, elemsToRefine, elemsToUnRefine, activeElements;

vector<double> refinementData, fluidProps, stagParams;
//vector<double> refinementData;
vector<double> fluidProps, stagParams;
vector<vector<int> > boundaryNodes, contElemData;
vector<vector<int> > NodeNumsAtLevel;
vector<vector<int> > DDconn;

vector<vector<double> > DirichletBCs, NeumannBCs, DerivativeBCs, pointBCs, Iconds;
vector<vector<double> > refineLimitVals, FluidOutputData;
vector<vector<double> > refineLimitVals, FluidOutputData, refinementData;

vector<node*> elems;
//typedef boost::ptr_vector<base> container;
Expand Down
58 changes: 50 additions & 8 deletions src/HBsplines/HBSplineBase1.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -314,7 +314,7 @@ void HBSplineBase::prepareInputData()
param.setZero();

// Refine the underlying grid if specified

/*
if(refinementData.size() > 0)
{
//printf("\n HBSplineBase : Refinement process ...... STARTED \n \n ");
Expand Down Expand Up @@ -348,12 +348,53 @@ void HBSplineBase::prepareInputData()
findUnique(elemsToRefine);
/*
printf(" elemsToRefine \n");
for(ii=0;ii<elemsToRefine.size();ii++)
printf("%5d \t", elemsToRefine[ii] );
printf("\n\n\n");
*/
MAX_LEVEL += 1;
applyRefinementProcess();
CURRENT_LEVEL += 1;
processBoundaryConditionsRefinedLevels();
//for(ii=0;ii<elemsToRefine.size();ii++)
//elems[elemsToRefine[ii]]->printSelf();
}
//printf("\n HBSplineBase : Refinement process ...... FINISHED \n \n ");
}
*/


if(refinementData.size() > 0)
{
//printf("\n HBSplineBase : Refinement process ...... STARTED \n \n ");

for(kk=0; kk<refinementData.size(); kk++)
{
for(int rr=0;rr<refinementData[kk][1];rr++)
{
elemsToRefine.clear();

switch( int (refinementData[kk][0]) )
{
default :
case 0:
refine(refinementData[kk][2]);
break;
case 1:
refinementforHemkerProblem();
break;
case 2:
refinementforAdvDiff1D();
break;
case 3:
refinementforAdvDiff2D();
break;
case 4:
pointBasedRefinement(refinementData[kk][2]);
break;
case 5:
limitBasedRefinement(rr);
break;
}

findUnique(elemsToRefine);

MAX_LEVEL += 1;
applyRefinementProcess();
Expand All @@ -365,6 +406,7 @@ void HBSplineBase::prepareInputData()
}
//printf("\n HBSplineBase : Refinement process ...... FINISHED \n \n ");
}
}

/////////////////////////////////
// assign boundary conditions to the element
Expand Down Expand Up @@ -620,7 +662,7 @@ void HBSplineBase::limitBasedRefinement(int kk)
double val[6];
myPoint knotBegin, knotIncr;

assert(refinementData[1] <= refineLimitVals.size());
//assert(refinementData[1] <= refineLimitVals.size());

if(DIM == 1)
{
Expand Down
16 changes: 9 additions & 7 deletions src/HBsplines/HBSplineCutFEM1.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,13 +213,16 @@ void HBSplineCutFEM::readInputData(std::ifstream &Ifile, MyString &line)

refinementData.resize(lvdTmp.n);

//cout << lvdTmp[0] << endl;
refinementData.resize(lvdTmp[0].n);
if(lvdTmp[0].n < 1)
prgError(2, fct, "invalid number of 'refinement type' !");
for(i=0;i<lvdTmp.n;i++)
{
if(lvdTmp[i].n < 1)
prgError(2, fct, "invalid number of 'refinement type' !");

for(j=0;j<lvdTmp[0].n;j++)
refinementData[j] = lvdTmp[0][j];
refinementData[i].resize(lvdTmp[i].n);

for(j=0;j<lvdTmp[i].n;j++)
refinementData[i][j] = lvdTmp[i][j];
}

//printVector(refinementData);

Expand All @@ -234,7 +237,6 @@ void HBSplineCutFEM::readInputData(std::ifstream &Ifile, MyString &line)

for(i=0;i<lvdTmp.n;i++)
{
//cout << lvdTmp[i] << endl;
refineLimitVals[i].resize(lvdTmp[i].n);
if(lvdTmp[i].n < 1)
prgError(2, fct, "invalid number of 'mesh refinement limits' !");
Expand Down
15 changes: 9 additions & 6 deletions src/HBsplines/HBSplineFEM1.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,13 +184,16 @@ void HBSplineFEM::readInputData(std::ifstream &Ifile, MyString &line)

refinementData.resize(lvdTmp.n);

//cout << lvdTmp[0] << endl;
refinementData.resize(lvdTmp[0].n);
if(lvdTmp[0].n < 1)
prgError(2, fct, "invalid number of 'refinement type' !");
for(i=0;i<lvdTmp.n;i++)
{
if(lvdTmp[1].n < 1)
prgError(2, fct, "invalid number of 'refinement type' !");

for(j=0;j<lvdTmp[0].n;j++)
refinementData[j] = lvdTmp[0][j];
refinementData[i].resize(lvdTmp[i].n);

for(j=0;j<lvdTmp[i].n;j++)
refinementData[i][j] = lvdTmp[i][j];
}

//printVector(refinementData);

Expand Down
33 changes: 32 additions & 1 deletion src/HBsplines/ImmersedFlexibleSolid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -412,6 +412,7 @@ void ImmersedFlexibleSolid::writeOutput()

int bb=0, type=0, nn2=0, dof=0, ii=0, kk=0;
double val_out=0.0;
VectorXd energy(9);

for(bb=0; bb<OutputData.size(); bb++)
{
Expand Down Expand Up @@ -470,9 +471,21 @@ void ImmersedFlexibleSolid::writeOutput()
sprintf(tmp," \t %12.6E", val_out);
break;

case 100 : // strain energy

val_out=0.0;
for(int ee=0; ee<nElem; ee++)
{
elems[ee]->computeEnergy(0, 0, energy);

val_out += energy[0];
}
sprintf(tmp," \t %12.6E", val_out);
break;

default :

prgError(1,"HBSplineFEM::writeImmersedSolidOutput()","invalid value of 'type'!");
prgError(1,"ImmersedFlexibleSolid::writeOutput()","invalid value of 'type'!");

break;
}
Expand Down Expand Up @@ -507,6 +520,7 @@ void ImmersedFlexibleSolid::postProcess(int index)
vtkSmartPointer<vtkFloatArray> vecVTK2 = vtkSmartPointer<vtkFloatArray>::New();
vtkSmartPointer<vtkFloatArray> vecVTK3 = vtkSmartPointer<vtkFloatArray>::New();
vtkSmartPointer<vtkFloatArray> vecVTK4 = vtkSmartPointer<vtkFloatArray>::New();
vtkSmartPointer<vtkFloatArray> cellVTK = vtkSmartPointer<vtkFloatArray>::New();

vtkSmartPointer<vtkXMLUnstructuredGridWriter> writerUGridVTK = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();

Expand All @@ -517,6 +531,7 @@ void ImmersedFlexibleSolid::postProcess(int index)
vecVTK2->SetName("velo");
vecVTK3->SetName("acce");
vecVTK4->SetName("force");
cellVTK->SetName("StrEngy");

vecVTK->SetNumberOfComponents(3);
vecVTK->SetNumberOfTuples(nNode);
Expand All @@ -527,6 +542,8 @@ void ImmersedFlexibleSolid::postProcess(int index)
vecVTK4->SetNumberOfComponents(3);
vecVTK4->SetNumberOfTuples(nNode);

cellVTK->SetNumberOfTuples(nElem);

vtkIdType pt[2];

if(DIM == 2)
Expand Down Expand Up @@ -671,11 +688,25 @@ void ImmersedFlexibleSolid::postProcess(int index)
}
}

// compute strain energy
VectorXd energy(9);
double val;
for(ee=0; ee<nElem; ee++)
{
elems[ee]->computeEnergy(0, 0, energy);

val = energy[0];

cellVTK->SetTuple1(ee, val);
}


uGridVTK->SetPoints(pointsVTK);
uGridVTK->GetPointData()->SetVectors(vecVTK);
uGridVTK->GetPointData()->AddArray(vecVTK2);
uGridVTK->GetPointData()->AddArray(vecVTK3);
uGridVTK->GetPointData()->AddArray(vecVTK4);
uGridVTK->GetCellData()->AddArray(cellVTK);

char fname[500];

Expand Down
77 changes: 47 additions & 30 deletions src/standardFEM/LagrangeElem2DBbarFbar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -699,7 +699,7 @@ void LagrangeElem2DBbarFbar::toPostprocess(int vartype, int varindex, int type,

void LagrangeElem2DBbarFbar::computeEnergy(int index1, int index2, VectorXd& energy)
{
// compute total energy for structural dynamic problems
// compute strain energy
///////////////////////////////////////////////////////////

elmDat = &(SolnData->ElemProp[elmType].data[0]);
Expand All @@ -711,16 +711,15 @@ void LagrangeElem2DBbarFbar::computeEnergy(int index1, int index2, VectorXd& en

VectorXd N(nlbf), dN_dx(nlbf), dN_dy(nlbf);

int err, isw, count, count1, index, ll = 0, ii, jj, gp1, gp2, TI, TIp1, TJ, TJp1;
int err, isw, count, count1, index, ll = 0, ii, jj, gp, TI, TIp1, TJ, TJp1;
int ind1, ind2, kk;

double rho = elmDat[5] ;

vector<double> gausspoints1, gaussweights1;

vector<double> gausspoints1, gausspoints2, gaussweights;
// compute determinant detF (tr(F)) in element centre

getGaussPoints1D(1, gausspoints1, gaussweights1);
getGaussPoints1D(1, gausspoints1, gaussweights);

param[0] = GeomData->gausspoints1[0];
param[1] = GeomData->gausspoints1[0];
Expand All @@ -742,28 +741,29 @@ void LagrangeElem2DBbarFbar::computeEnergy(int index1, int index2, VectorXd& en

//cout << " finite = " << finite << endl;

int nGP = (int) elmDat[0] ;

getGaussPoints1D(nGP, gausspoints1, gaussweights1);
int nGP = 4 ;
getGaussPointsQuad(nGP, gausspoints1, gausspoints2, gaussweights);

energy.setZero();
energy.setZero();

for(gp2=0;gp2<nGP;gp2++)
{
JacMult = gaussweights1[gp2] * thick;
MatrixXd FF(3,3), Finv(3,3), sigma(3,3), LagTens(3,3), Iden(3,3), PK2(3,3);
FF.setZero();
PK2.setZero();
sigma.setZero();
Iden.setZero();
Iden(0,0) = 1.0; Iden(1,1) = 1.0; Iden(2,2) = 1.0;

param[1] = gausspoints1[gp2];

for(gp1=0;gp1<nGP;gp1++)
{
param[0] = gausspoints1[gp1];
for(gp=0; gp<nGP; gp++)
{
param[0] = gausspoints1[gp];
param[1] = gausspoints2[gp];

GeomData->computeBasisFunctions2D(0, 2, degree, param, nodeNums, &N(0), &dN_dx(0), &dN_dy(0), Jac);

fact = gaussweights1[gp1] * JacMult;

fact = gaussweights[gp] * thick;
dvol0 = Jac * fact;
dvol = dvol0;
dvol = dvol0;

//if(axsy)
//dvol *= 2.0*PI*yy;
Expand Down Expand Up @@ -822,21 +822,38 @@ void LagrangeElem2DBbarFbar::computeEnergy(int index1, int index2, VectorXd& en
strain[2] = F33 - 1.0;
strain[3] = 0.5 * (F[1] + F[2]);

//printf(" %14.12f \t %14.12f \t %14.12f \t %14.12f \n ", strain[0], strain[1], strain[2], strain[3]);
//printf(" %14.12f \t %14.12f \t %14.12f \t %14.12f \n ", stress[0], stress[1], stress[2], stress[3]);
//printf(" dvol = %14.12f \n\n ", dvol);

velo[0] = computeValueDot(0, N);
velo[1] = computeValueDot(1, N);
//velo[0] = computeValueDot(0, N);
//velo[1] = computeValueDot(1, N);

// kinetic energy
energy[0] += (0.5*dvol0*rho*(velo[0]*velo[0]+velo[1]*velo[1]));
energy[1] += (0.5*dvol*(strain[0]*stress[0]+strain[1]*stress[1]+strain[2]*stress[2]+strain[3]*stress[3]));
//energy[0] += (0.5*dvol0*rho*(velo[0]*velo[0]+velo[1]*velo[1]));
//energy[1] += (0.5*dvol*(strain[0]*stress[0]+strain[1]*stress[1]+strain[2]*stress[2]+strain[3]*stress[3]));

FF(0,0) = F[0];
FF(0,1) = F[2];
FF(1,0) = F[1];
FF(1,1) = F[3];
FF(2,2) = 1.0;
Finv = FF.inverse();

LagTens = 0.5*(FF.transpose()*FF - Iden);

sigma(0,0) = stress[0];
sigma(1,1) = stress[1];
sigma(2,2) = stress[2];

sigma(0,1) = stress[3]; sigma(1,0) = stress[3];

}//gp1
}//gp2
PK2 = detF*Finv*sigma*Finv.transpose();

//printf("\n\n");
// strain energy
energy[0] += (0.5*dvol* ( PK2(0,0)*LagTens(0,0) + PK2(1,1)*LagTens(1,1) + PK2(2,2)*LagTens(2,2) + 2.0*PK2(0,1)*LagTens(0,1) ) );

//printf(" %14.12f \t %14.12f \t %14.12f \t %14.12f \n ", strain[0], strain[1], strain[2], strain[3]);
//printf(" %14.12f \t %14.12f \t %14.12f \t %14.12f \n ", stress[0], stress[1], stress[2], stress[3]);
//printf(" dvol = %14.12f \n\n ", dvol);
}//gp
//printf("\n\n");

return;
}
Expand Down
Loading

0 comments on commit c52d701

Please sign in to comment.