Skip to content

Commit

Permalink
Storing DS now
Browse files Browse the repository at this point in the history
  • Loading branch information
@axla-io
axla-io committed Jul 3, 2023
1 parent 710c52e commit ade25bd
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 16 deletions.
7 changes: 4 additions & 3 deletions examples/prom/nonlinear_elasticity_global_rom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,6 +164,7 @@ class RomOperator : public TimeDependentOperator
std::vector<double> eqp_rw;
std::vector<int> eqp_qp;
Vector eqp_coef;
Vector eqp_DS_coef;
const bool fastIntegration = true;

int rank;
Expand Down Expand Up @@ -1735,7 +1736,7 @@ RomOperator::RomOperator(HyperelasticOperator *fom_,
<< fom->fespace.GetNE() << endl;

GetEQPCoefficients_HyperelasticNLFIntegrator(&(fom->fespace), eqp_rw, eqp_qp,
ir_eqp, model, V_v, rank, eqp_coef);
ir_eqp, model, V_v, rank, eqp_coef, eqp_DS_coef);
}
}

Expand Down Expand Up @@ -1771,7 +1772,7 @@ void RomOperator::Mult_Hyperreduced(const Vector &vx, Vector &dvx_dt) const
HyperelasticNLFIntegrator_ComputeReducedEQP_Fast(&(fom->fespace),eqp_rw,
eqp_qp, ir_eqp, model,
x0, V_x, V_v, x_librom,
eqp_coef, rank, resEQP);
eqp_coef, eqp_DS_coef, rank, resEQP);
}
else
HyperelasticNLFIntegrator_ComputeReducedEQP(&(fom->fespace), eqp_rw,
Expand Down Expand Up @@ -1909,5 +1910,5 @@ void RomOperator::SetEQP(CAROM::Vector *eqpSol)
<< fom->fespace.GetNE() << endl;

GetEQPCoefficients_HyperelasticNLFIntegrator(&(fom->fespace), eqp_rw, eqp_qp,
ir_eqp, model, V_v, rank, eqp_coef);
ir_eqp, model, V_v, rank, eqp_coef, eqp_DS_coef);
}
48 changes: 35 additions & 13 deletions examples/prom/nonlinear_elasticity_global_rom_eqp.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ class HyperelasticOperator : public TimeDependentOperator
void GetEQPCoefficients_HyperelasticNLFIntegrator(ParFiniteElementSpace *fesR,
std::vector<double> const &rw, std::vector<int> const &qp,
const IntegrationRule *ir, NeoHookeanModel *model,
CAROM::Matrix const &V_v, const int rank, Vector &coef)
CAROM::Matrix const &V_v, const int rank, Vector &coef, Vector &DS_coef)
{
const int rvdim = V_v.numColumns();
const int fomdim = V_v.numRows();
Expand Down Expand Up @@ -96,6 +96,17 @@ void GetEQPCoefficients_HyperelasticNLFIntegrator(ParFiniteElementSpace *fesR,
coef.SetSize(elvect_size * rw.size() * rvdim);
coef = 0.0;

// Vector for storing DS
const FiniteElement *fe = fesR->GetFE(0);
dof = fe->GetDof();
dim = fe->GetDim();
DenseMatrix DSh(dof, dim);
DenseMatrix DS(dof, dim);
DenseMatrix Jrt(dim);
DS_coef.SetSize(dof * dim * rw.size() * rvdim);
DS_coef = 0.0;
int index = 0;

// For every basis vector
for (int j = 0; j < rvdim; ++j)
{
Expand Down Expand Up @@ -139,6 +150,19 @@ void GetEQPCoefficients_HyperelasticNLFIntegrator(ParFiniteElementSpace *fesR,
{
coef[k + (i * elvect_size) + (j * rw.size() * elvect_size)] = vj_e[k] * rw[i] * t;
}

// Calculate DS and store
CalcInverse(eltrans->Jacobian(), Jrt);
fe->CalcDShape(ip, DSh);
Mult(DSh, Jrt, DS);
for (int ii = 0; ii < dof; ++ii)
{
for (int jj = 0; jj < dim; ++jj)
{
index = jj + ii * dim;
DS_coef[index + (i * dof*dim) + (j * rw.size() * dof*dim)] = DS.Elem(ii,jj);
}
}
}
}
}
Expand Down Expand Up @@ -282,7 +306,7 @@ void HyperelasticNLFIntegrator_ComputeReducedEQP(ParFiniteElementSpace *fesR,
void HyperelasticNLFIntegrator_ComputeReducedEQP_Fast(ParFiniteElementSpace *fesR,
std::vector<double> const &rw, std::vector<int> const &qp, const IntegrationRule *ir, NeoHookeanModel *model,
const Vector *x0, CAROM::Matrix const &V_x, CAROM::Matrix const &V_v, CAROM::Vector const &x,
Vector const &coef, const int rank, Vector &res)
Vector const &coef, Vector const &DS_coef, const int rank, Vector &res)
{

const int rxdim = V_x.numColumns();
Expand Down Expand Up @@ -332,9 +356,8 @@ void HyperelasticNLFIntegrator_ComputeReducedEQP_Fast(ParFiniteElementSpace *fes
fe = fesR->GetFE(0);
dof = fe->GetDof();
dim = fe->GetDim();
DenseMatrix DSh(dof, dim);
int index = 0;
DenseMatrix DS(dof, dim);
DenseMatrix Jrt(dim);
DenseMatrix Jpt(dim);
DenseMatrix P_f(dim);
DenseMatrix PMatI; // Extract element dofs
Expand Down Expand Up @@ -384,18 +407,17 @@ void HyperelasticNLFIntegrator_ComputeReducedEQP_Fast(ParFiniteElementSpace *fes
eltrans->SetIntPoint(&ip);
model->SetTransformation(*eltrans);

// Get the transformation weight
// double t = eltrans->Weight();
// t = 1.0; //Hack

// Compute action of nonlinear operator
CalcInverse(eltrans->Jacobian(), Jrt); // TODO: incorporate in coefficient calculation, if it makes sense
fe->CalcDShape(ip, DSh); // TODO: incorporate in coefficient calculation, if it makes sense
Mult(DSh, Jrt, DS); // TODO: incorporate in coefficient calculation, if it makes sense
for (int ii = 0; ii < dof; ++ii)
{
for (int jj = 0; jj < dim; ++jj)
{
index = jj + ii * dim;
DS.Elem(ii,jj) = DS_coef[index + (i * elvect.Size()) + (j * qp.size() * elvect.Size())];
}
}

MultAtB(PMatI, DS, Jpt);
model->EvalP(Jpt, P_f);
// P_f *= (t * rw[i]); // NB: Not by ip.weight
AddMultABt(DS, P_f, PMatO);

// Calculate r[i] = ve_j^T * elvect
Expand Down

0 comments on commit ade25bd

Please sign in to comment.