52#ifdef NEKTAR_USE_PETSC
75 int n = geom->GetNumVerts(), i, j;
81 for (i = 0; i < n; ++i)
83 for (j = 0; j < n - 1; ++j)
85 map(j, i) = (*geom->GetVertex(i))[j];
120 for (i = 0; i < n - 1; ++i)
122 for (j = 0; j < n - 1; ++j)
124 mapred(i, j) = newmap(i, j);
152 const int nVel =
m_fields[0]->GetCoordim(0);
155 ASSERTL0(nVel > 1,
"Linear elastic solver not set up for"
156 " this dimension (only 2D/3D supported).");
166 std::dynamic_pointer_cast<MultiRegions::ContField>(
m_fields[0]);
173 const int nEl =
m_fields[0]->GetExpSize();
186 for (n = 0; n < nEl; ++n)
189 sizeBnd[n] = nVel * exp->NumBndryCoeffs();
190 sizeInt[n] = nVel * exp->GetNcoeffs() - sizeBnd[n];
191 exp->GetBoundaryMap(
m_bmap[n]);
192 exp->GetInteriorMap(
m_imap[n]);
198 sizeBnd, sizeBnd, blkmatStorage);
230 const int nEl =
m_fields[0]->GetExpSize();
231 const int nVel =
m_fields[0]->GetCoordim(0);
243 bool verbose =
m_session->DefinesCmdLineArgument(
"verbose");
244 bool root =
m_session->GetComm()->GetRank() == 0;
249 for (n = 0; n < nEl; ++n)
252 const int nPhys = exp->GetTotPoints();
263 exp->DetShapeType(), *exp,
factors,
266 exp->DetShapeType(), *exp,
factors,
274 mat[0][0] = exp->GenMatrix(matkeyA);
276 mat[1][0] = mat[0][1];
277 mat[1][1] = exp->GenMatrix(matkeyD);
284 std::cout <<
"\rBuilding matrix system: "
285 << (int)(100.0 * n / nEl) <<
"%" << std::flush;
291 for (n = 0; n < nEl; ++n)
294 const int nPhys = exp->GetTotPoints();
310 exp->DetShapeType(), *exp,
factors,
313 exp->DetShapeType(), *exp,
factors,
316 exp->DetShapeType(), *exp,
factors,
325 mat[0][0] = exp->GenMatrix(matkeyA);
329 mat[1][0] = mat[0][1];
330 mat[1][1] = exp->GenMatrix(matkeyE);
333 mat[2][0] = mat[0][2];
334 mat[2][1] = mat[1][2];
335 mat[2][2] = exp->GenMatrix(matkeyI);
342 std::cout <<
"\rBuilding matrix system: "
343 << (int)(100.0 * n / nEl) <<
"%" << std::flush;
350 std::cout <<
"\rBuilding matrix system: done." << std::endl;
382 const int nVel =
m_fields[0]->GetCoordim(0);
411#ifdef NEKTAR_USE_PETSC
432 const int nCoeffs =
m_fields[0]->GetNcoeffs();
443 std::string tempEval;
444 m_session->LoadSolverInfo(
"Temperature", tempEval,
"None");
446 if (tempEval ==
"Jacobian")
451 for (nv = 0; nv < nVel; ++nv)
457 for (i = 0; i <
m_fields[0]->GetExpSize(); ++i)
464 int offset =
m_fields[0]->GetPhys_Offset(i);
466 if (exp->GetMetricInfo()->GetGtype() ==
481 else if (tempEval ==
"Metric")
484 m_graph->GetAllQuadGeoms().size() == 0) ||
486 m_graph->GetAllPrismGeoms().size() == 0 &&
487 m_graph->GetAllPyrGeoms().size() == 0 &&
488 m_graph->GetAllHexGeoms().size() == 0),
489 "LinearIdealMetric temperature only implemented for "
490 "two-dimensional triangular meshes or three-dimensional "
491 "tetrahedral meshes.");
496 for (nv = 0; nv < nVel; ++nv)
502 for (i = 0; i < nVel; ++i)
509 for (i = 0; i <
m_fields[0]->GetExpSize(); ++i)
515 exp->GetMetricInfo()->GetDeriv(pkey);
516 int offset =
m_fields[0]->GetPhys_Offset(i);
525 for (j = 0; j < deriv[0][0].size(); ++j)
527 for (k = 0; k < nVel; ++k)
529 for (l = 0; l < nVel; ++l)
531 jac(l, k) = deriv[k][l][j];
535 jacIdeal = jac * i2rm;
540 char jobvl =
'N', jobvr =
'V';
541 int worklen = 8 * nVel, info;
551 &(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
552 &(evec.GetPtr())[0], nVel, &work[0], worklen,
559 for (nv = 0; nv < nVel; ++nv)
561 eval(nv, nv) =
m_beta * (
sqrt(eval(nv, nv)) - 1.0);
566 for (k = 0; k < nVel; ++k)
568 for (l = 0; l < nVel; ++l)
575 if (deriv[0][0].size() != exp->GetTotPoints())
578 for (k = 0; k < nVel; ++k)
580 for (l = 0; l < nVel; ++l)
591 for (i = 0; i < nVel; ++i)
593 for (j = 0; j < nVel; ++j)
604 else if (tempEval !=
"None")
606 ASSERTL0(
false,
"Unknown temperature form: " + tempEval);
622 for (nv = 0; nv < nVel; ++nv)
625 m_fields[nv]->IProductWRTBase(forcing[nv], tmp);
629 m_fields[nv]->ImposeDirichletConditions(loc_inout);
633 for (i = 0; i <
m_fields[nv]->GetExpSize(); ++i)
635 int nBnd =
m_bmap[i].size();
636 int nInt =
m_imap[i].size();
637 int offset =
m_fields[nv]->GetCoeff_Offset(i);
639 for (j = 0; j < nBnd; ++j)
641 forCoeffs[nVel * offset + nv * nBnd + j] =
642 -1.0 * tmp[offset +
m_bmap[i][j]];
643 inout[nVel * offset + nv * nBnd + j] =
644 loc_inout[offset +
m_bmap[i][j]];
646 for (j = 0; j < nInt; ++j)
648 forCoeffs[nVel * (offset + nBnd) + nv * nInt + j] =
649 -1.0 * tmp[offset +
m_imap[i][j]];
665 for (nv = 0; nv < nVel; ++nv)
667 for (i = 0; i <
m_fields[nv]->GetExpSize(); ++i)
669 int nBnd =
m_bmap[i].size();
670 int nInt =
m_imap[i].size();
671 int offset =
m_fields[nv]->GetCoeff_Offset(i);
673 for (j = 0; j < nBnd; ++j)
676 inout[nVel * offset + nv * nBnd + j];
678 for (j = 0; j < nInt; ++j)
681 inout[nVel * (offset + nBnd) + nv * nInt + j];
706 const int nVel = mat.GetRows();
707 const int nB = exp->NumBndryCoeffs();
708 const int nI = exp->GetNcoeffs() - nB;
709 const int nBnd = exp->NumBndryCoeffs() * nVel;
710 const int nInt = exp->GetNcoeffs() * nVel - nBnd;
723 for (i = 0; i < nVel; ++i)
725 for (j = 0; j < nVel; ++j)
728 for (k = 0; k < nB; ++k)
730 for (l = 0; l < nB; ++l)
732 (*A)(k + i * nB, l + j * nB) =
736 for (l = 0; l < nI; ++l)
738 (*B)(k + i * nB, l + j * nI) =
744 for (k = 0; k < nI; ++k)
746 for (l = 0; l < nB; ++l)
748 (*C)(k + i * nI, l + j * nB) =
752 for (l = 0; l < nI; ++l)
754 (*D)(k + i * nI, l + j * nI) =
764 (*A) = (*A) - (*B) * (*C);
785 const int k1,
const int k2,
const NekDouble scale,
788 const int nCoeffs = exp->GetNcoeffs();
789 const int nPhys = exp->GetTotPoints();
798 for (i = 0; i < nCoeffs; ++i)
803 exp->BwdTrans(tmp1, tmp2);
804 exp->PhysDeriv(k1, tmp2, tmp3);
805 exp->IProductWRTDerivBase(k2, tmp3, tmp1);
808 &(ret->GetPtr())[0] + i * nCoeffs, 1);
816 std::vector<std::string> &variables)
818 const int nVel =
m_fields[0]->GetCoordim(0);
819 const int nCoeffs =
m_fields[0]->GetNcoeffs();
826 for (
int i = 0; i < nVel; ++i)
830 fieldcoeffs.push_back(tFwd);
831 variables.push_back(
"ThermStressDiv" + std::to_string(i));
839 for (
int i = 0; i < nVel; ++i)
841 for (
int j = 0; j < nVel; ++j)
845 fieldcoeffs.push_back(tFwd);
846 variables.push_back(
"ThermStress" + std::to_string(i) +
#define ASSERTL0(condition, msg)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Array< OneD, Array< OneD, unsigned int > > m_imap
Interior maps for each of the fields.
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
DNekScalBlkMatSharedPtr m_Dinv
Matrix of elemental components.
NekDouble m_nu
Poisson ratio.
void v_InitObject(bool DeclareFields=true) override
Set up the linear elasticity system.
void BuildMatrixSystem()
Build matrix system for linear elasticity equations.
LinearElasticSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Default constructor.
void SetStaticCondBlock(const int n, const LocalRegions::ExpansionSharedPtr exp, Array< TwoD, DNekMatSharedPtr > &mat)
Given a block matrix for an element, construct its static condensation matrices.
Array< OneD, Array< OneD, NekDouble > > m_temperature
Storage for the temperature terms.
CoupledAssemblyMapSharedPtr m_assemblyMap
Assembly map for the coupled (u,v,w) system.
void v_DoSolve() override
Solve elliptic linear elastic system.
DNekScalBlkMatSharedPtr m_BinvD
Matrix of elemental components.
NekDouble m_beta
Parameter dictating amount of thermal stress to add.
DNekScalBlkMatSharedPtr m_schurCompl
Schur complement boundary-boundary matrix.
void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables) override
NekDouble m_E
Young's modulus.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_stress
Storage for the thermal stress terms.
static std::string className
Name of class.
Array< OneD, Array< OneD, unsigned int > > m_bmap
Boundary maps for each of the fields.
DNekMatSharedPtr BuildLaplacianIJMatrix(const int k1, const int k2, const NekDouble scale, LocalRegions::ExpansionSharedPtr exp)
Construct a LaplacianIJ matrix for a given expansion.
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Generate summary at runtime.
DNekScalBlkMatSharedPtr m_C
Matrix of elemental components.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Describe a linear system.
A base class for describing how to solve specific equations.
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareFeld=true)
Initialisation object for EquationSystem.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
static void Dgeev(const char &uplo, const char &lrev, const int &n, const double *a, const int &lda, double *wr, double *wi, double *rev, const int &ldr, double *lev, const int &ldv, double *work, const int &lwork, int &info)
Solve general real matrix eigenproblem.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
@ beta
Gauss Radau pinned at x=-1,.
std::shared_ptr< Expansion > ExpansionSharedPtr
static PreconditionerSharedPtr NullPreconditionerSharedPtr
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
std::shared_ptr< ContField > ContFieldSharedPtr
std::vector< std::pair< std::string, std::string > > SummaryList
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< Geometry > GeometrySharedPtr
@ eLinearAdvectionReaction
std::map< ConstFactorType, NekDouble > ConstFactorMap
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
StdRegions::ConstFactorMap factors
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekMat > DNekMatSharedPtr
DNekMat MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom)
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)