Nektar++
Static Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | Friends | List of all members
Nektar::LinearElasticSystem Class Reference

Base class for linear elastic system. More...

#include <LinearElasticSystem.h>

Inheritance diagram for Nektar::LinearElasticSystem:
[legend]

Static Public Member Functions

static EquationSystemSharedPtr create (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Creates an instance of this class. More...
 

Static Public Attributes

static std::string className
 Name of class. More...
 

Protected Member Functions

 LinearElasticSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Default constructor. More...
 
void BuildMatrixSystem ()
 Build matrix system for linear elasticity equations. More...
 
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. More...
 
DNekMatSharedPtr BuildLaplacianIJMatrix (const int k1, const int k2, const NekDouble scale, LocalRegions::ExpansionSharedPtr exp)
 Construct a LaplacianIJ matrix for a given expansion. More...
 
void v_InitObject (bool DeclareFields=true) override
 Set up the linear elasticity system. More...
 
void v_GenerateSummary (SolverUtils::SummaryList &s) override
 Generate summary at runtime. More...
 
void v_DoSolve () override
 Solve elliptic linear elastic system. More...
 
void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables) override
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members. More...
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareFeld=true)
 Initialisation object for EquationSystem. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true)
 Virtual function for initialisation implementation. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoSolve ()
 Virtual function for solve implementation. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Virtual function for the L_inf error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the L_2 error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_H1Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the H_1 error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys ()
 Virtual function for transformation to physical space. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff ()
 Virtual function for transformation to coefficient space. More...
 
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &l)
 Virtual function for generating summary information. More...
 
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution (unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 
virtual SOLVER_UTILS_EXPORT void v_Output (void)
 
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure (void)
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp (void)
 Virtual function to identify if operator is negated in DoSolve. More...
 
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 

Protected Attributes

NekDouble m_nu
 Poisson ratio. More...
 
NekDouble m_E
 Young's modulus. More...
 
NekDouble m_beta
 Parameter dictating amount of thermal stress to add. More...
 
CoupledAssemblyMapSharedPtr m_assemblyMap
 Assembly map for the coupled (u,v,w) system. More...
 
DNekScalBlkMatSharedPtr m_schurCompl
 Schur complement boundary-boundary matrix. More...
 
DNekScalBlkMatSharedPtr m_BinvD
 Matrix of elemental \( B^{-1}D \) components. More...
 
DNekScalBlkMatSharedPtr m_C
 Matrix of elemental \( C \) components. More...
 
DNekScalBlkMatSharedPtr m_Dinv
 Matrix of elemental \( D^{-1} \) components. More...
 
Array< OneD, Array< OneD, unsigned int > > m_bmap
 Boundary maps for each of the fields. More...
 
Array< OneD, Array< OneD, unsigned int > > m_imap
 Interior maps for each of the fields. More...
 
Array< OneD, Array< OneD, NekDouble > > m_temperature
 Storage for the temperature terms. More...
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_stress
 Storage for the thermal stress terms. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
bool m_verbose
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtrm_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_fields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
NekDouble m_lastCheckTime
 
NekDouble m_TimeIncrementFactor
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_infosteps
 Number of time steps between outputting status information. More...
 
int m_iterPIT = 0
 Number of parallel-in-time time iteration. More...
 
int m_windowPIT = 0
 Index of windows for parallel-in-time time iteration. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
Array< OneD, NekDoublem_movingFrameData
 Moving reference frame status in the inertial frame X, Y, Z, Theta_x, Theta_y, Theta_z, U, V, W, Omega_x, Omega_y, Omega_z, A_x, A_y, A_z, DOmega_x, DOmega_y, DOmega_z, pivot_x, pivot_y, pivot_z. More...
 
std::vector< std::string > m_strFrameData
 variable name in m_movingFrameData More...
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error. More...
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous) More...
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous) More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
int m_HomoDirec
 number of homogenous directions More...
 

Friends

class MemoryManager< LinearElasticSystem >
 Class may only be instantiated through the MemoryManager. More...
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void InitObject (bool DeclareField=true)
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise (bool dumpInitialConditions=true)
 Perform any initialisation necessary before solving the problem. More...
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem. More...
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space. More...
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space. More...
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve. More...
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name. More...
 
template<class T >
std::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name. More...
 
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession ()
 Get Session name. More...
 
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure ()
 Get pressure field if available. More...
 
SOLVER_UTILS_EXPORT void ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 
SOLVER_UTILS_EXPORT void PrintSummary (std::ostream &out)
 Print a summary of parameters and solver characteristics. More...
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda. More...
 
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction (std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
 Get a SessionFunction by name. More...
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields. More...
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the L2 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields. More...
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT NekDouble H1Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the H1 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDoubleErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write checkpoint file of custom data fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write input fields to the given filename. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
 Input field data from the given file to multiple domains. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, std::vector< std::string > &fieldStr, Array< OneD, Array< OneD, NekDouble > > &coeffs)
 Output a field. Input field data into array from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, MultiRegions::ExpListSharedPtr &pField, std::string &pFieldName)
 Output a field. Input field data into ExpList from the given file. More...
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary. More...
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMapUpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, const Array< OneD, const NekDouble > &input)
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > & UpdatePhysField (const int i)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT int GetInfoSteps ()
 
SOLVER_UTILS_EXPORT void SetInfoSteps (int num)
 
SOLVER_UTILS_EXPORT void SetIterationNumberPIT (int num)
 
SOLVER_UTILS_EXPORT void SetWindowNumberPIT (int num)
 
SOLVER_UTILS_EXPORT Array< OneD, const Array< OneD, NekDouble > > GetTraceNormals ()
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
SOLVER_UTILS_EXPORT bool NegatedOp ()
 Identify if operator is negated in DoSolve. More...
 
- Protected Types inherited from Nektar::SolverUtils::EquationSystem
enum  HomogeneousType { eHomogeneous1D , eHomogeneous2D , eHomogeneous3D , eNotHomogeneous }
 Parameter for homogeneous expansions. More...
 
- Static Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
static std::string equationSystemTypeLookupIds []
 
static std::string projectionTypeLookupIds []
 

Detailed Description

Base class for linear elastic system.

Definition at line 49 of file LinearElasticSystem.h.

Constructor & Destructor Documentation

◆ LinearElasticSystem()

Nektar::LinearElasticSystem::LinearElasticSystem ( const LibUtilities::SessionReaderSharedPtr pSession,
const SpatialDomains::MeshGraphSharedPtr pGraph 
)
protected

Default constructor.

Definition at line 134 of file LinearElasticSystem.cpp.

137 : EquationSystem(pSession, pGraph)
138{
139}
SOLVER_UTILS_EXPORT EquationSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises EquationSystem class members.

Member Function Documentation

◆ BuildLaplacianIJMatrix()

DNekMatSharedPtr Nektar::LinearElasticSystem::BuildLaplacianIJMatrix ( const int  k1,
const int  k2,
const NekDouble  scale,
LocalRegions::ExpansionSharedPtr  exp 
)
protected

Construct a LaplacianIJ matrix for a given expansion.

This routine constructs matrices whose entries contain the evaluation of

\[ \partial_{\tt k1} \phi_i \partial_{\tt k2} \phi_j \,dx \]

Definition at line 784 of file LinearElasticSystem.cpp.

787{
788 const int nCoeffs = exp->GetNcoeffs();
789 const int nPhys = exp->GetTotPoints();
790 int i;
791
792 DNekMatSharedPtr ret =
793 MemoryManager<DNekMat>::AllocateSharedPtr(nCoeffs, nCoeffs, 0.0, eFULL);
794
795 Array<OneD, NekDouble> tmp2(nPhys);
796 Array<OneD, NekDouble> tmp3(nPhys);
797
798 for (i = 0; i < nCoeffs; ++i)
799 {
800 Array<OneD, NekDouble> tmp1(nCoeffs, 0.0);
801 tmp1[i] = 1.0;
802
803 exp->BwdTrans(tmp1, tmp2);
804 exp->PhysDeriv(k1, tmp2, tmp3);
805 exp->IProductWRTDerivBase(k2, tmp3, tmp1);
806
807 Vmath::Smul(nCoeffs, scale, &tmp1[0], 1,
808 &(ret->GetPtr())[0] + i * nCoeffs, 1);
809 }
810
811 return ret;
812}
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, and Vmath::Smul().

Referenced by BuildMatrixSystem().

◆ BuildMatrixSystem()

void Nektar::LinearElasticSystem::BuildMatrixSystem ( )
protected

Build matrix system for linear elasticity equations.

This routine constructs the matrix discretisation arising from the weak formulation of the linear elasticity equations. The resulting matrix system is then passed to LinearElasticSystem::SetStaticCondBlock in order to construct the statically condensed system.

All of the matrices involved in the construction of the divergence of the stress tensor are Laplacian-like. We use the variable coefficient functionality within the library to multiply the Laplacian by the appropriate constants for all diagonal terms, and off-diagonal terms use LinearElasticSystem::BuildLaplacianIJMatrix to construct matrices containing

\[ \int \partial_{x_i} \phi_k \partial_{x_j} \phi_l \]

where mixed derivative terms are present. Symmetry (in terms of k,l) is exploited to avoid constructing this matrix repeatedly.

Todo:
Make static condensation optional and construct full system instead.

Definition at line 228 of file LinearElasticSystem.cpp.

229{
230 const int nEl = m_fields[0]->GetExpSize();
231 const int nVel = m_fields[0]->GetCoordim(0);
232
234 int n;
235
236 // Factors map for matrix keys.
238
239 // Calculate various constants
240 NekDouble mu = m_E * 0.5 / (1.0 + m_nu);
241 NekDouble lambda = m_E * m_nu / (1.0 + m_nu) / (1.0 - 2.0 * m_nu);
242
243 bool verbose = m_session->DefinesCmdLineArgument("verbose");
244 bool root = m_session->GetComm()->GetRank() == 0;
245
246 // Loop over each element and construct matrices.
247 if (nVel == 2)
248 {
249 for (n = 0; n < nEl; ++n)
250 {
251 exp = m_fields[0]->GetExp(n);
252 const int nPhys = exp->GetTotPoints();
253 Array<OneD, NekDouble> aArr(nPhys, lambda + 2 * mu);
254 Array<OneD, NekDouble> bArr(nPhys, mu);
255
256 StdRegions::VarCoeffMap varcoeffA, varcoeffD;
257 varcoeffA[StdRegions::eVarCoeffD00] = aArr;
258 varcoeffA[StdRegions::eVarCoeffD11] = bArr;
259 varcoeffD[StdRegions::eVarCoeffD00] = bArr;
260 varcoeffD[StdRegions::eVarCoeffD11] = aArr;
261
262 LocalRegions::MatrixKey matkeyA(StdRegions::eLaplacian,
263 exp->DetShapeType(), *exp, factors,
264 varcoeffA);
265 LocalRegions::MatrixKey matkeyD(StdRegions::eLaplacian,
266 exp->DetShapeType(), *exp, factors,
267 varcoeffD);
268
269 /*
270 * mat holds the linear operator [ A B ] acting on [ u ].
271 * [ C D ] [ v ]
272 */
273 Array<TwoD, DNekMatSharedPtr> mat(2, 2);
274 mat[0][0] = exp->GenMatrix(matkeyA);
275 mat[0][1] = BuildLaplacianIJMatrix(1, 0, mu + lambda, exp);
276 mat[1][0] = mat[0][1];
277 mat[1][1] = exp->GenMatrix(matkeyD);
278
279 // Set up the statically condensed block for this element.
280 SetStaticCondBlock(n, exp, mat);
281
282 if (verbose && root)
283 {
284 std::cout << "\rBuilding matrix system: "
285 << (int)(100.0 * n / nEl) << "%" << std::flush;
286 }
287 }
288 }
289 else if (nVel == 3)
290 {
291 for (n = 0; n < nEl; ++n)
292 {
293 exp = m_fields[0]->GetExp(n);
294 const int nPhys = exp->GetTotPoints();
295 Array<OneD, NekDouble> aArr(nPhys, lambda + 2 * mu);
296 Array<OneD, NekDouble> bArr(nPhys, mu);
297
298 StdRegions::VarCoeffMap varcoeffA, varcoeffE, varcoeffI;
299 varcoeffA[StdRegions::eVarCoeffD00] = aArr;
300 varcoeffA[StdRegions::eVarCoeffD11] = bArr;
301 varcoeffA[StdRegions::eVarCoeffD22] = bArr;
302 varcoeffE[StdRegions::eVarCoeffD00] = bArr;
303 varcoeffE[StdRegions::eVarCoeffD11] = aArr;
304 varcoeffE[StdRegions::eVarCoeffD22] = bArr;
305 varcoeffI[StdRegions::eVarCoeffD00] = bArr;
306 varcoeffI[StdRegions::eVarCoeffD11] = bArr;
307 varcoeffI[StdRegions::eVarCoeffD22] = aArr;
308
309 LocalRegions::MatrixKey matkeyA(StdRegions::eLaplacian,
310 exp->DetShapeType(), *exp, factors,
311 varcoeffA);
312 LocalRegions::MatrixKey matkeyE(StdRegions::eLaplacian,
313 exp->DetShapeType(), *exp, factors,
314 varcoeffE);
315 LocalRegions::MatrixKey matkeyI(StdRegions::eLaplacian,
316 exp->DetShapeType(), *exp, factors,
317 varcoeffI);
318
319 /*
320 * mat holds the linear operator [ A B C ] acting on [ u ].
321 * [ D E F ] [ v ]
322 * [ G H I ] [ w ]
323 */
324 Array<TwoD, DNekMatSharedPtr> mat(3, 3);
325 mat[0][0] = exp->GenMatrix(matkeyA);
326 mat[0][1] = BuildLaplacianIJMatrix(1, 0, mu + lambda, exp);
327 mat[0][2] = BuildLaplacianIJMatrix(2, 0, mu + lambda, exp);
328
329 mat[1][0] = mat[0][1];
330 mat[1][1] = exp->GenMatrix(matkeyE);
331 mat[1][2] = BuildLaplacianIJMatrix(2, 1, mu + lambda, exp);
332
333 mat[2][0] = mat[0][2];
334 mat[2][1] = mat[1][2];
335 mat[2][2] = exp->GenMatrix(matkeyI);
336
337 // Set up the statically condensed block for this element.
338 SetStaticCondBlock(n, exp, mat);
339
340 if (verbose && root)
341 {
342 std::cout << "\rBuilding matrix system: "
343 << (int)(100.0 * n / nEl) << "%" << std::flush;
344 }
345 }
346 }
347
348 if (verbose && root)
349 {
350 std::cout << "\rBuilding matrix system: done." << std::endl;
351 }
352}
NekDouble m_nu
Poisson ratio.
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.
NekDouble m_E
Young's modulus.
DNekMatSharedPtr BuildLaplacianIJMatrix(const int k1, const int k2, const NekDouble scale, LocalRegions::ExpansionSharedPtr exp)
Construct a LaplacianIJ matrix for a given expansion.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition: StdRegions.hpp:375
StdRegions::ConstFactorMap factors
double NekDouble

References BuildLaplacianIJMatrix(), Nektar::StdRegions::eLaplacian, Nektar::StdRegions::eVarCoeffD00, Nektar::StdRegions::eVarCoeffD11, Nektar::StdRegions::eVarCoeffD22, Nektar::VarcoeffHashingTest::factors, m_E, Nektar::SolverUtils::EquationSystem::m_fields, m_nu, Nektar::SolverUtils::EquationSystem::m_session, and SetStaticCondBlock().

Referenced by v_DoSolve().

◆ create()

static EquationSystemSharedPtr Nektar::LinearElasticSystem::create ( const LibUtilities::SessionReaderSharedPtr pSession,
const SpatialDomains::MeshGraphSharedPtr pGraph 
)
inlinestatic

Creates an instance of this class.

Definition at line 56 of file LinearElasticSystem.h.

59 {
62 pGraph);
63 p->InitObject();
64 return p;
65 }
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

◆ SetStaticCondBlock()

void Nektar::LinearElasticSystem::SetStaticCondBlock ( const int  n,
const LocalRegions::ExpansionSharedPtr  exp,
Array< TwoD, DNekMatSharedPtr > &  mat 
)
protected

Given a block matrix for an element, construct its static condensation matrices.

This routine essentially duplicates the logic present in many of the LocalRegions matrix generation routines to construct the statically condensed equivalent of mat to pass to the GlobalLinSys solver.

Parameters
nElement/block number.
expPointer to expansion.
matBlock matrix containing matrix operator.

Definition at line 701 of file LinearElasticSystem.cpp.

704{
705 int i, j, k, l;
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;
711 const MatrixStorage s = eFULL; // Maybe look into doing symmetric
712 // version of this?
713
722
723 for (i = 0; i < nVel; ++i)
724 {
725 for (j = 0; j < nVel; ++j)
726 {
727 // Boundary-boundary and boundary-interior
728 for (k = 0; k < nB; ++k)
729 {
730 for (l = 0; l < nB; ++l)
731 {
732 (*A)(k + i * nB, l + j * nB) =
733 (*mat[i][j])(m_bmap[n][k], m_bmap[n][l]);
734 }
735
736 for (l = 0; l < nI; ++l)
737 {
738 (*B)(k + i * nB, l + j * nI) =
739 (*mat[i][j])(m_bmap[n][k], m_imap[n][l]);
740 }
741 }
742
743 // Interior-boundary / interior-interior
744 for (k = 0; k < nI; ++k)
745 {
746 for (l = 0; l < nB; ++l)
747 {
748 (*C)(k + i * nI, l + j * nB) =
749 (*mat[i][j])(m_imap[n][k], m_bmap[n][l]);
750 }
751
752 for (l = 0; l < nI; ++l)
753 {
754 (*D)(k + i * nI, l + j * nI) =
755 (*mat[i][j])(m_imap[n][k], m_imap[n][l]);
756 }
757 }
758 }
759 }
760
761 // Construct static condensation matrices.
762 D->Invert();
763 (*B) = (*B) * (*D);
764 (*A) = (*A) - (*B) * (*C);
765
766 DNekScalMatSharedPtr tmp_mat;
767 m_schurCompl->SetBlock(
769 m_BinvD->SetBlock(
770 n, n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, B));
771 m_C->SetBlock(
772 n, n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, C));
773 m_Dinv->SetBlock(
774 n, n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, D));
775}
Array< OneD, Array< OneD, unsigned int > > m_imap
Interior maps for each of the fields.
DNekScalBlkMatSharedPtr m_Dinv
Matrix of elemental components.
DNekScalBlkMatSharedPtr m_BinvD
Matrix of elemental components.
DNekScalBlkMatSharedPtr m_schurCompl
Schur complement boundary-boundary matrix.
Array< OneD, Array< OneD, unsigned int > > m_bmap
Boundary maps for each of the fields.
DNekScalBlkMatSharedPtr m_C
Matrix of elemental components.
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, m_BinvD, m_bmap, m_C, m_Dinv, m_imap, and m_schurCompl.

Referenced by BuildMatrixSystem().

◆ v_DoSolve()

void Nektar::LinearElasticSystem::v_DoSolve ( void  )
overrideprotectedvirtual

Solve elliptic linear elastic system.

The solve proceeds as follows:

  • Create a MultiRegions::GlobalLinSys object.
  • Evaluate a forcing function from the session file.
  • If the temperature term is enabled, evaluate this and add to forcing.
  • Apply Dirichlet boundary conditions.
  • Scatter forcing into the correct ordering according to the coupled assembly map.
  • Do the solve.
  • Scatter solution back to fields and backwards transform to physical space.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 379 of file LinearElasticSystem.cpp.

380{
381 int i, j, k, l, nv;
382 const int nVel = m_fields[0]->GetCoordim(0);
383
384 // Build initial matrix system.
386
387 // Now we've got the matrix system set up, create a GlobalLinSys
388 // object. We mask ourselves as LinearAdvectionReaction to create a full
389 // matrix instead of symmetric storage.
390 MultiRegions::GlobalLinSysKey key(StdRegions::eLinearAdvectionReaction,
393
394 // Currently either direct or iterative static condensation is
395 // supported.
396 if (m_assemblyMap->GetGlobalSysSolnType() ==
398 {
402 }
403 else if (m_assemblyMap->GetGlobalSysSolnType() ==
405 {
410 }
411#ifdef NEKTAR_USE_PETSC
412 else if (m_assemblyMap->GetGlobalSysSolnType() ==
414 {
418 }
419#endif
420#ifdef NEKTAR_USE_MPI
421 else if (m_assemblyMap->GetGlobalSysSolnType() ==
423 {
427 }
428#endif
429
430 linSys->Initialise(m_assemblyMap);
431
432 const int nCoeffs = m_fields[0]->GetNcoeffs();
433
434 //
435 // -- Evaluate forcing functions
436 //
437
438 // Evaluate the forcing function from the XML file.
439 Array<OneD, Array<OneD, NekDouble>> forcing(nVel);
440 GetFunction("Forcing")->Evaluate(forcing);
441
442 // Add temperature term
443 std::string tempEval;
444 m_session->LoadSolverInfo("Temperature", tempEval, "None");
445
446 if (tempEval == "Jacobian")
447 {
448 // Allocate storage
449 m_temperature = Array<OneD, Array<OneD, NekDouble>>(nVel);
450
451 for (nv = 0; nv < nVel; ++nv)
452 {
453 Array<OneD, NekDouble> tmp;
454 m_temperature[nv] =
455 Array<OneD, NekDouble>(m_fields[nv]->GetNpoints());
456
457 for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
458 {
459 // Calculate element area
460 LocalRegions::ExpansionSharedPtr exp = m_fields[0]->GetExp(i);
461 LibUtilities::PointsKeyVector pkey = exp->GetPointsKeys();
462 Array<OneD, NekDouble> jac = exp->GetMetricInfo()->GetJac(pkey);
463
464 int offset = m_fields[0]->GetPhys_Offset(i);
465
466 if (exp->GetMetricInfo()->GetGtype() ==
468 {
469 Vmath::Smul(exp->GetTotPoints(), m_beta, jac, 1,
470 tmp = m_temperature[nv] + offset, 1);
471 }
472 else
473 {
474 Vmath::Fill(exp->GetTotPoints(), m_beta * jac[0],
475 tmp = m_temperature[nv] + offset, 1);
476 }
477 }
478 m_fields[nv]->PhysDeriv(nv, m_temperature[nv], forcing[nv]);
479 }
480 }
481 else if (tempEval == "Metric")
482 {
483 ASSERTL0((m_fields[0]->GetCoordim(0) == 2 &&
484 m_graph->GetAllQuadGeoms().size() == 0) ||
485 (m_fields[0]->GetCoordim(0) == 3 &&
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.");
492
493 m_temperature = Array<OneD, Array<OneD, NekDouble>>(nVel);
494 m_stress = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nVel);
495
496 for (nv = 0; nv < nVel; ++nv)
497 {
498 m_temperature[nv] =
499 Array<OneD, NekDouble>(m_fields[nv]->GetNpoints(), 0.0);
500 m_stress[nv] = Array<OneD, Array<OneD, NekDouble>>(nVel);
501
502 for (i = 0; i < nVel; ++i)
503 {
504 m_stress[nv][i] =
505 Array<OneD, NekDouble>(m_fields[nv]->GetNpoints(), 0.0);
506 }
507 }
508
509 for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
510 {
511 // Calculate element area
512 LocalRegions::ExpansionSharedPtr exp = m_fields[0]->GetExp(i);
513 LibUtilities::PointsKeyVector pkey = exp->GetPointsKeys();
514 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> deriv =
515 exp->GetMetricInfo()->GetDeriv(pkey);
516 int offset = m_fields[0]->GetPhys_Offset(i);
517
518 DNekMat i2rm = MappingIdealToRef(exp->GetGeom());
519
520 // Compute metric tensor
521 DNekMat jac(nVel, nVel, 0.0, eFULL);
522 DNekMat jacIdeal(nVel, nVel, 0.0, eFULL);
523 DNekMat metric(nVel, nVel, 0.0, eFULL);
524
525 for (j = 0; j < deriv[0][0].size(); ++j)
526 {
527 for (k = 0; k < nVel; ++k)
528 {
529 for (l = 0; l < nVel; ++l)
530 {
531 jac(l, k) = deriv[k][l][j];
532 }
533 }
534
535 jacIdeal = jac * i2rm;
536 metric = Transpose(jacIdeal) * jacIdeal;
537
538 // Compute eigenvalues/eigenvectors of metric tensor using
539 // ideal mapping.
540 char jobvl = 'N', jobvr = 'V';
541 int worklen = 8 * nVel, info;
542
543 DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
544 DNekMat evec(nVel, nVel, 0.0, eFULL);
545 DNekMat evecinv(nVel, nVel, 0.0, eFULL);
546 Array<OneD, NekDouble> vl(nVel * nVel);
547 Array<OneD, NekDouble> work(worklen);
548 Array<OneD, NekDouble> wi(nVel);
549
550 Lapack::Dgeev(jobvl, jobvr, nVel, metric.GetRawPtr(), nVel,
551 &(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
552 &(evec.GetPtr())[0], nVel, &work[0], worklen,
553 info);
554
555 evecinv = evec;
556 evecinv.Invert();
557
558 // rescaling of the eigenvalues
559 for (nv = 0; nv < nVel; ++nv)
560 {
561 eval(nv, nv) = m_beta * (sqrt(eval(nv, nv)) - 1.0);
562 }
563
564 DNekMat beta = evec * eval * evecinv;
565
566 for (k = 0; k < nVel; ++k)
567 {
568 for (l = 0; l < nVel; ++l)
569 {
570 m_stress[k][l][offset + j] = beta(k, l);
571 }
572 }
573 }
574
575 if (deriv[0][0].size() != exp->GetTotPoints())
576 {
577 Array<OneD, NekDouble> tmp;
578 for (k = 0; k < nVel; ++k)
579 {
580 for (l = 0; l < nVel; ++l)
581 {
582 Vmath::Fill(exp->GetTotPoints(), m_stress[k][l][offset],
583 tmp = m_stress[k][l] + offset, 1);
584 }
585 }
586 }
587 }
588
589 // Calculate divergence of stress tensor.
590 Array<OneD, NekDouble> tmpderiv(m_fields[0]->GetNpoints());
591 for (i = 0; i < nVel; ++i)
592 {
593 for (j = 0; j < nVel; ++j)
594 {
595 m_fields[i]->PhysDeriv(j, m_stress[i][j], tmpderiv);
596 Vmath::Vadd(m_fields[i]->GetNpoints(), tmpderiv, 1,
597 m_temperature[i], 1, m_temperature[i], 1);
598 }
599
601 forcing[i], 1);
602 }
603 }
604 else if (tempEval != "None")
605 {
606 ASSERTL0(false, "Unknown temperature form: " + tempEval);
607 }
608
609 // Set up some temporary storage.
610 //
611 // - forCoeffs holds the forcing coefficients in a local ordering;
612 // however note that the ordering is different and dictated by the
613 // assembly map. We loop over each element, then the boundary degrees of
614 // freedom for u, boundary for v, followed by the interior for u and then
615 // interior for v.
616 // - inout holds the Dirichlet degrees of freedom in the local ordering,
617 // which have the boundary conditions imposed
618 Array<OneD, NekDouble> forCoeffs(nVel * nCoeffs, 0.0);
619 Array<OneD, NekDouble> inout(nVel * nCoeffs, 0.0);
620 Array<OneD, NekDouble> tmp(nCoeffs);
621
622 for (nv = 0; nv < nVel; ++nv)
623 {
624 // Take the inner product of the forcing function.
625 m_fields[nv]->IProductWRTBase(forcing[nv], tmp);
626
627 // Impose Dirichlet condition on field which should be initialised
628 Array<OneD, NekDouble> loc_inout = m_fields[nv]->UpdateCoeffs();
629 m_fields[nv]->ImposeDirichletConditions(loc_inout);
630
631 // Scatter forcing into RHS vector according to the ordering dictated in
632 // the comment above.
633 for (i = 0; i < m_fields[nv]->GetExpSize(); ++i)
634 {
635 int nBnd = m_bmap[i].size();
636 int nInt = m_imap[i].size();
637 int offset = m_fields[nv]->GetCoeff_Offset(i);
638
639 for (j = 0; j < nBnd; ++j)
640 {
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]];
645 }
646 for (j = 0; j < nInt; ++j)
647 {
648 forCoeffs[nVel * (offset + nBnd) + nv * nInt + j] =
649 -1.0 * tmp[offset + m_imap[i][j]];
650 }
651 }
652 }
653
654 //
655 // -- Perform solve
656 //
657 // Solve.
658 linSys->Solve(forCoeffs, inout, m_assemblyMap);
659
660 //
661 // -- Postprocess
662 //
663
664 // Scatter back to field degrees of freedom
665 for (nv = 0; nv < nVel; ++nv)
666 {
667 for (i = 0; i < m_fields[nv]->GetExpSize(); ++i)
668 {
669 int nBnd = m_bmap[i].size();
670 int nInt = m_imap[i].size();
671 int offset = m_fields[nv]->GetCoeff_Offset(i);
672
673 for (j = 0; j < nBnd; ++j)
674 {
675 m_fields[nv]->UpdateCoeffs()[offset + m_bmap[i][j]] =
676 inout[nVel * offset + nv * nBnd + j];
677 }
678 for (j = 0; j < nInt; ++j)
679 {
680 m_fields[nv]->UpdateCoeffs()[offset + m_imap[i][j]] =
681 inout[nVel * (offset + nBnd) + nv * nInt + j];
682 }
683 }
684 m_fields[nv]->BwdTrans(m_fields[nv]->GetCoeffs(),
685 m_fields[nv]->UpdatePhys());
686 }
687}
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
void BuildMatrixSystem()
Build matrix system for linear elasticity equations.
Array< OneD, Array< OneD, NekDouble > > m_temperature
Storage for the temperature terms.
CoupledAssemblyMapSharedPtr m_assemblyMap
Assembly map for the coupled (u,v,w) system.
NekDouble m_beta
Parameter dictating amount of thermal stress to add.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_stress
Storage for the thermal stress terms.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
SOLVER_UTILS_EXPORT int GetNpoints()
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.
Definition: Lapack.hpp:329
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
static PreconditionerSharedPtr NullPreconditionerSharedPtr
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:51
@ eDeformed
Geometry is curved or has non-constant factors.
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:50
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
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.
Definition: Vmath.hpp:180
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::beta, BuildMatrixSystem(), Lapack::Dgeev(), Nektar::SpatialDomains::eDeformed, Nektar::eDIAGONAL, Nektar::MultiRegions::eDirectStaticCond, Nektar::eFULL, Nektar::MultiRegions::eIterativeStaticCond, Nektar::StdRegions::eLinearAdvectionReaction, Nektar::MultiRegions::ePETScStaticCond, Nektar::MultiRegions::eXxtStaticCond, Vmath::Fill(), Nektar::SolverUtils::EquationSystem::GetFunction(), Nektar::SolverUtils::EquationSystem::GetNpoints(), m_assemblyMap, m_beta, m_BinvD, m_bmap, m_C, m_Dinv, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_graph, m_imap, m_schurCompl, Nektar::SolverUtils::EquationSystem::m_session, m_stress, m_temperature, Nektar::MappingIdealToRef(), Nektar::MultiRegions::NullPreconditionerSharedPtr, Vmath::Smul(), tinysimd::sqrt(), Nektar::Transpose(), Vmath::Vadd(), and Vmath::Vcopy().

Referenced by Nektar::IterativeElasticSystem::v_DoSolve().

◆ v_ExtraFldOutput()

void Nektar::LinearElasticSystem::v_ExtraFldOutput ( std::vector< Array< OneD, NekDouble > > &  fieldcoeffs,
std::vector< std::string > &  variables 
)
overrideprotectedvirtual

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 814 of file LinearElasticSystem.cpp.

817{
818 const int nVel = m_fields[0]->GetCoordim(0);
819 const int nCoeffs = m_fields[0]->GetNcoeffs();
820
821 if (m_temperature.size() == 0)
822 {
823 return;
824 }
825
826 for (int i = 0; i < nVel; ++i)
827 {
828 Array<OneD, NekDouble> tFwd(nCoeffs);
829 m_fields[i]->FwdTrans(m_temperature[i], tFwd);
830 fieldcoeffs.push_back(tFwd);
831 variables.push_back("ThermStressDiv" + std::to_string(i));
832 }
833
834 if (m_stress.size() == 0)
835 {
836 return;
837 }
838
839 for (int i = 0; i < nVel; ++i)
840 {
841 for (int j = 0; j < nVel; ++j)
842 {
843 Array<OneD, NekDouble> tFwd(nCoeffs);
844 m_fields[i]->FwdTrans(m_stress[i][j], tFwd);
845 fieldcoeffs.push_back(tFwd);
846 variables.push_back("ThermStress" + std::to_string(i) +
847 std::to_string(j));
848 }
849 }
850}

References Nektar::SolverUtils::EquationSystem::m_fields, m_stress, and m_temperature.

◆ v_GenerateSummary()

void Nektar::LinearElasticSystem::v_GenerateSummary ( SolverUtils::SummaryList s)
overrideprotectedvirtual

Generate summary at runtime.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 357 of file LinearElasticSystem.cpp.

358{
360
361 AddSummaryItem(s, "Young's modulus", m_E);
362 AddSummaryItem(s, "Poisson ratio", m_nu);
363}
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47

References Nektar::SolverUtils::AddSummaryItem(), m_E, m_nu, and Nektar::SolverUtils::EquationSystem::SessionSummary().

Referenced by Nektar::IterativeElasticSystem::v_GenerateSummary().

◆ v_InitObject()

void Nektar::LinearElasticSystem::v_InitObject ( bool  DeclareFields = true)
overrideprotectedvirtual

Set up the linear elasticity system.

This routine loads the E and nu variables from the session file, creates a coupled assembly map and creates containers for the statically condensed block matrix system.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 148 of file LinearElasticSystem.cpp.

149{
150 EquationSystem::v_InitObject(DeclareFields);
151
152 const int nVel = m_fields[0]->GetCoordim(0);
153 int n;
154
155 ASSERTL0(nVel > 1, "Linear elastic solver not set up for"
156 " this dimension (only 2D/3D supported).");
157
158 // Make sure that we have Young's modulus and Poisson ratio set.
159 m_session->LoadParameter("E", m_E, 1.00);
160 m_session->LoadParameter("nu", m_nu, 0.25);
161 m_session->LoadParameter("Beta", m_beta, 1.00);
162
163 // Create a coupled assembly map which allows us to tie u, v and w
164 // fields together.
166 std::dynamic_pointer_cast<MultiRegions::ContField>(m_fields[0]);
168 m_session, m_graph, u->GetLocalToGlobalMap(),
169 m_fields[0]->GetBndConditions(), m_fields);
170
171 // Figure out size of our new matrix systems by looping over all expansions
172 // and multiply number of coefficients by velocity components.
173 const int nEl = m_fields[0]->GetExpSize();
175
176 Array<OneD, unsigned int> sizeBnd(nEl);
177 Array<OneD, unsigned int> sizeInt(nEl);
178
179 // Allocate storage for boundary and interior map caches.
180 m_bmap = Array<OneD, Array<OneD, unsigned int>>(nEl);
181 m_imap = Array<OneD, Array<OneD, unsigned int>>(nEl);
182
183 // Cache interior and boundary maps to avoid recalculating
184 // constantly. Should really be handled by manager in StdRegions but not
185 // really working yet.
186 for (n = 0; n < nEl; ++n)
187 {
188 exp = m_fields[0]->GetExp(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]);
193 }
194
195 // Create block matrix storage for the statically condensed system.
196 MatrixStorage blkmatStorage = eDIAGONAL;
198 sizeBnd, sizeBnd, blkmatStorage);
200 blkmatStorage);
202 blkmatStorage);
204 blkmatStorage);
205}
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareFeld=true)
Initialisation object for EquationSystem.
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::eDIAGONAL, m_assemblyMap, m_beta, m_BinvD, m_bmap, m_C, m_Dinv, m_E, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_graph, m_imap, m_nu, m_schurCompl, Nektar::SolverUtils::EquationSystem::m_session, and Nektar::SolverUtils::EquationSystem::v_InitObject().

Referenced by Nektar::IterativeElasticSystem::v_InitObject().

Friends And Related Function Documentation

◆ MemoryManager< LinearElasticSystem >

friend class MemoryManager< LinearElasticSystem >
friend

Class may only be instantiated through the MemoryManager.

Definition at line 1 of file LinearElasticSystem.h.

Member Data Documentation

◆ className

std::string Nektar::LinearElasticSystem::className
static
Initial value:
=
"LinearElasticSystem", LinearElasticSystem::create)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
EquationSystemFactory & GetEquationSystemFactory()

Name of class.

Definition at line 68 of file LinearElasticSystem.h.

◆ m_assemblyMap

CoupledAssemblyMapSharedPtr Nektar::LinearElasticSystem::m_assemblyMap
protected

Assembly map for the coupled (u,v,w) system.

Definition at line 78 of file LinearElasticSystem.h.

Referenced by v_DoSolve(), and v_InitObject().

◆ m_beta

NekDouble Nektar::LinearElasticSystem::m_beta
protected

Parameter dictating amount of thermal stress to add.

Definition at line 76 of file LinearElasticSystem.h.

Referenced by v_DoSolve(), and v_InitObject().

◆ m_BinvD

DNekScalBlkMatSharedPtr Nektar::LinearElasticSystem::m_BinvD
protected

Matrix of elemental \( B^{-1}D \) components.

Definition at line 82 of file LinearElasticSystem.h.

Referenced by SetStaticCondBlock(), v_DoSolve(), and v_InitObject().

◆ m_bmap

Array<OneD, Array<OneD, unsigned int> > Nektar::LinearElasticSystem::m_bmap
protected

Boundary maps for each of the fields.

Definition at line 88 of file LinearElasticSystem.h.

Referenced by SetStaticCondBlock(), v_DoSolve(), and v_InitObject().

◆ m_C

DNekScalBlkMatSharedPtr Nektar::LinearElasticSystem::m_C
protected

Matrix of elemental \( C \) components.

Definition at line 84 of file LinearElasticSystem.h.

Referenced by SetStaticCondBlock(), v_DoSolve(), and v_InitObject().

◆ m_Dinv

DNekScalBlkMatSharedPtr Nektar::LinearElasticSystem::m_Dinv
protected

Matrix of elemental \( D^{-1} \) components.

Definition at line 86 of file LinearElasticSystem.h.

Referenced by SetStaticCondBlock(), v_DoSolve(), and v_InitObject().

◆ m_E

NekDouble Nektar::LinearElasticSystem::m_E
protected

Young's modulus.

Definition at line 74 of file LinearElasticSystem.h.

Referenced by BuildMatrixSystem(), v_GenerateSummary(), and v_InitObject().

◆ m_imap

Array<OneD, Array<OneD, unsigned int> > Nektar::LinearElasticSystem::m_imap
protected

Interior maps for each of the fields.

Definition at line 90 of file LinearElasticSystem.h.

Referenced by SetStaticCondBlock(), v_DoSolve(), and v_InitObject().

◆ m_nu

NekDouble Nektar::LinearElasticSystem::m_nu
protected

Poisson ratio.

Definition at line 72 of file LinearElasticSystem.h.

Referenced by BuildMatrixSystem(), v_GenerateSummary(), and v_InitObject().

◆ m_schurCompl

DNekScalBlkMatSharedPtr Nektar::LinearElasticSystem::m_schurCompl
protected

Schur complement boundary-boundary matrix.

Definition at line 80 of file LinearElasticSystem.h.

Referenced by SetStaticCondBlock(), v_DoSolve(), and v_InitObject().

◆ m_stress

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::LinearElasticSystem::m_stress
protected

Storage for the thermal stress terms.

Definition at line 94 of file LinearElasticSystem.h.

Referenced by v_DoSolve(), and v_ExtraFldOutput().

◆ m_temperature

Array<OneD, Array<OneD, NekDouble> > Nektar::LinearElasticSystem::m_temperature
protected

Storage for the temperature terms.

Definition at line 92 of file LinearElasticSystem.h.

Referenced by v_DoSolve(), and v_ExtraFldOutput().