Nektar++
Public Member Functions | 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]

Public Member Functions

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...
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void SetUpTraceNormals (void)
 
SOLVER_UTILS_EXPORT void InitObject ()
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise ()
 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 NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. 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 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 GetFinalTime ()
 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 NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
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 void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp ()
 Virtual function to identify if operator is negated in DoSolve. More...
 

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...
 
virtual void v_InitObject ()
 Set up the linear elasticity system. More...
 
virtual void v_GenerateSummary (SolverUtils::SummaryList &s)
 Generate summary at runtime. More...
 
virtual void v_DoSolve ()
 Solve elliptic linear elastic system. More...
 
virtual void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 
- 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_DoInitialise ()
 Virtual function for initialisation 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 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_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)
 

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
 
bool m_root
 
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_timestepMax = -1.0
 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_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...
 
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

- 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 []
 

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 126 of file LinearElasticSystem.cpp.

129  : EquationSystem(pSession, pGraph)
130 {
131 }
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 
)

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 791 of file LinearElasticSystem.cpp.

796 {
797  const int nCoeffs = exp->GetNcoeffs();
798  const int nPhys = exp->GetTotPoints();
799  int i;
800 
802  nCoeffs, nCoeffs, 0.0, eFULL);
803 
804  Array<OneD, NekDouble> tmp2(nPhys);
805  Array<OneD, NekDouble> tmp3(nPhys);
806 
807  for (i = 0; i < nCoeffs; ++i)
808  {
809  Array<OneD, NekDouble> tmp1(nCoeffs, 0.0);
810  tmp1[i] = 1.0;
811 
812  exp->BwdTrans ( tmp1, tmp2);
813  exp->PhysDeriv (k1, tmp2, tmp3);
814  exp->IProductWRTDerivBase(k2, tmp3, tmp1);
815 
816  Vmath::Smul(
817  nCoeffs, scale, &tmp1[0], 1, &(ret->GetPtr())[0]+i*nCoeffs, 1);
818  }
819 
820  return ret;
821 }
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:69
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.cpp:225

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

Referenced by BuildMatrixSystem().

◆ BuildMatrixSystem()

void Nektar::LinearElasticSystem::BuildMatrixSystem ( )

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 223 of file LinearElasticSystem.cpp.

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

References BuildLaplacianIJMatrix(), Nektar::StdRegions::eLaplacian, Nektar::StdRegions::eVarCoeffD00, Nektar::StdRegions::eVarCoeffD11, Nektar::StdRegions::eVarCoeffD22, 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  {
60  EquationSystemSharedPtr p = MemoryManager<
61  LinearElasticSystem>::AllocateSharedPtr(pSession, pGraph);
62  p->InitObject();
63  return p;
64  }
LinearElasticSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Default constructor.
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.

References CellMLToNektar.cellml_metadata::p.

◆ SetStaticCondBlock()

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

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 703 of file LinearElasticSystem.cpp.

707 {
708  int i, j, k, l;
709  const int nVel = mat.GetRows();
710  const int nB = exp->NumBndryCoeffs();
711  const int nI = exp->GetNcoeffs() - nB;
712  const int nBnd = exp->NumBndryCoeffs() * nVel;
713  const int nInt = exp->GetNcoeffs() * nVel - nBnd;
714  const MatrixStorage s = eFULL; // Maybe look into doing symmetric
715  // version of this?
716 
718  MemoryManager<DNekMat>::AllocateSharedPtr(nBnd, nBnd, 0.0, s);
719  DNekMatSharedPtr B =
720  MemoryManager<DNekMat>::AllocateSharedPtr(nBnd, nInt, 0.0, s);
721  DNekMatSharedPtr C =
722  MemoryManager<DNekMat>::AllocateSharedPtr(nInt, nBnd, 0.0, s);
723  DNekMatSharedPtr D =
724  MemoryManager<DNekMat>::AllocateSharedPtr(nInt, nInt, 0.0, s);
725 
726  for (i = 0; i < nVel; ++i)
727  {
728  for (j = 0; j < nVel; ++j)
729  {
730  // Boundary-boundary and boundary-interior
731  for (k = 0; k < nB; ++k)
732  {
733  for (l = 0; l < nB; ++l)
734  {
735  (*A)(k + i*nB, l + j*nB) =
736  (*mat[i][j])(m_bmap[n][k], m_bmap[n][l]);
737  }
738 
739  for (l = 0; l < nI; ++l)
740  {
741  (*B)(k + i*nB, l + j*nI) =
742  (*mat[i][j])(m_bmap[n][k], m_imap[n][l]);
743  }
744  }
745 
746  // Interior-boundary / interior-interior
747  for (k = 0; k < nI; ++k)
748  {
749  for (l = 0; l < nB; ++l)
750  {
751  (*C)(k + i*nI, l + j*nB) =
752  (*mat[i][j])(m_imap[n][k], m_bmap[n][l]);
753  }
754 
755  for (l = 0; l < nI; ++l)
756  {
757  (*D)(k + i*nI, l + j*nI) =
758  (*mat[i][j])(m_imap[n][k], m_imap[n][l]);
759  }
760  }
761  }
762  }
763 
764  // Construct static condensation matrices.
765  D->Invert();
766  (*B) = (*B)*(*D);
767  (*A) = (*A) - (*B)*(*C);
768 
769  DNekScalMatSharedPtr tmp_mat;
770  m_schurCompl->SetBlock(
772  1.0, A));
773  m_BinvD ->SetBlock(
775  1.0, B));
776  m_C ->SetBlock(
778  1.0, C));
779  m_Dinv ->SetBlock(
781  1.0, D));
782 }
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  )
protectedvirtual

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.

Reimplemented in Nektar::IterativeElasticSystem.

Definition at line 374 of file LinearElasticSystem.cpp.

375 {
376  int i, j, k, l, nv;
377  const int nVel = m_fields[0]->GetCoordim(0);
378 
379  // Build initial matrix system.
381 
382  // Now we've got the matrix system set up, create a GlobalLinSys
383  // object. We mask ourselves as LinearAdvectionReaction to create a full
384  // matrix instead of symmetric storage.
385  MultiRegions::GlobalLinSysKey key(
388 
389  // Currently either direct or iterative static condensation is
390  // supported.
391  if (m_assemblyMap->GetGlobalSysSolnType() == MultiRegions::eDirectStaticCond)
392  {
393  linSys = MemoryManager<
394  MultiRegions::GlobalLinSysDirectStaticCond>::AllocateSharedPtr(
395  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
396  m_assemblyMap);
397  }
398  else if (m_assemblyMap->GetGlobalSysSolnType() ==
400  {
401  linSys = MemoryManager<
402  MultiRegions::GlobalLinSysIterativeStaticCond>::AllocateSharedPtr(
403  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
405  }
406 #ifdef NEKTAR_USE_PETSC
407  else if (m_assemblyMap->GetGlobalSysSolnType() ==
409  {
410  linSys = MemoryManager<
411  MultiRegions::GlobalLinSysPETScStaticCond>::AllocateSharedPtr(
412  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
413  m_assemblyMap);
414  }
415 #endif
416 #ifdef NEKTAR_USE_MPI
417  else if (m_assemblyMap->GetGlobalSysSolnType() ==
419  {
420  linSys = MemoryManager<
421  MultiRegions::GlobalLinSysXxtStaticCond>::AllocateSharedPtr(
422  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
423  m_assemblyMap);
424  }
425 #endif
426 
427  linSys->Initialise(m_assemblyMap);
428 
429  const int nCoeffs = m_fields[0]->GetNcoeffs();
430 
431  //
432  // -- Evaluate forcing functions
433  //
434 
435  // Evaluate the forcing function from the XML file.
436  Array<OneD, Array<OneD, NekDouble> > forcing(nVel);
437  GetFunction("Forcing")->Evaluate(forcing);
438 
439  // Add temperature term
440  string tempEval;
441  m_session->LoadSolverInfo("Temperature", tempEval, "None");
442 
443  if (tempEval == "Jacobian")
444  {
445  // Allocate storage
446  m_temperature = Array<OneD, Array<OneD, NekDouble> >(nVel);
447 
448  for (nv = 0; nv < nVel; ++nv)
449  {
450  Array<OneD, NekDouble> tmp;
451  m_temperature[nv] = Array<OneD, NekDouble>(
452  m_fields[nv]->GetNpoints());
453 
454  for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
455  {
456  // Calculate element area
458  m_fields[0]->GetExp(i);
460  exp->GetPointsKeys();
461  Array<OneD, NekDouble> jac =
462  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] = Array<OneD, NekDouble>(
499  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] = Array<OneD, NekDouble>(
505  m_fields[nv]->GetNpoints(), 0.0);
506  }
507  }
508 
509  for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
510  {
511  // Calculate element area
513  m_fields[0]->GetExp(i);
514  LibUtilities::PointsKeyVector pkey = exp->GetPointsKeys();
515  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > deriv =
516  exp->GetMetricInfo()->GetDeriv(pkey);
517  int offset = m_fields[0]->GetPhys_Offset(i);
518 
519  DNekMat i2rm = MappingIdealToRef(exp->GetGeom());
520 
521  // Compute metric tensor
522  DNekMat jac (nVel, nVel, 0.0, eFULL);
523  DNekMat jacIdeal(nVel, nVel, 0.0, eFULL);
524  DNekMat metric (nVel, nVel, 0.0, eFULL);
525 
526  for (j = 0; j < deriv[0][0].size(); ++j)
527  {
528  for (k = 0; k < nVel; ++k)
529  {
530  for (l = 0; l < nVel; ++l)
531  {
532  jac(l,k) = deriv[k][l][j];
533  }
534  }
535 
536  jacIdeal = jac * i2rm;
537  metric = Transpose(jacIdeal) * jacIdeal;
538 
539  // Compute eigenvalues/eigenvectors of metric tensor using
540  // ideal mapping.
541  char jobvl = 'N', jobvr = 'V';
542  int worklen = 8*nVel, info;
543 
544  DNekMat eval (nVel, nVel, 0.0, eDIAGONAL);
545  DNekMat evec (nVel, nVel, 0.0, eFULL);
546  DNekMat evecinv(nVel, nVel, 0.0, eFULL);
547  Array<OneD, NekDouble> vl (nVel*nVel);
548  Array<OneD, NekDouble> work(worklen);
549  Array<OneD, NekDouble> wi (nVel);
550 
551  Lapack::Dgeev(jobvl, jobvr, nVel, metric.GetRawPtr(), nVel,
552  &(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
553  &(evec.GetPtr())[0], nVel,
554  &work[0], worklen, info);
555 
556  evecinv = evec;
557  evecinv.Invert();
558 
559  // rescaling of the eigenvalues
560  for (nv = 0; nv < nVel; ++nv)
561  {
562  eval(nv,nv) = m_beta * (sqrt(eval(nv,nv)) - 1.0);
563  }
564 
565  DNekMat beta = evec * eval * evecinv;
566 
567  for (k = 0; k < nVel; ++k)
568  {
569  for (l = 0; l < nVel; ++l)
570  {
571  m_stress[k][l][offset+j] = beta(k,l);
572  }
573  }
574  }
575 
576  if (deriv[0][0].size() != exp->GetTotPoints())
577  {
578  Array<OneD, NekDouble> tmp;
579  for (k = 0; k < nVel; ++k)
580  {
581  for (l = 0; l < nVel; ++l)
582  {
583  Vmath::Fill(
584  exp->GetTotPoints(), m_stress[k][l][offset],
585  tmp = m_stress[k][l] + offset, 1);
586  }
587  }
588  }
589  }
590 
591  // Calculate divergence of stress tensor.
592  Array<OneD, NekDouble> tmpderiv(m_fields[0]->GetNpoints());
593  for (i = 0; i < nVel; ++i)
594  {
595  for (j = 0; j < nVel; ++j)
596  {
597  m_fields[i]->PhysDeriv(j, m_stress[i][j], tmpderiv);
598  Vmath::Vadd (m_fields[i]->GetNpoints(), tmpderiv, 1,
599  m_temperature[i], 1, m_temperature[i], 1);
600  }
601 
603  m_temperature[i], 1, forcing[i], 1);
604  }
605  }
606  else if (tempEval != "None")
607  {
608  ASSERTL0(false, "Unknown temperature form: " + tempEval);
609  }
610 
611  // Set up some temporary storage.
612  //
613  // - forCoeffs holds the forcing coefficients in a local ordering;
614  // however note that the ordering is different and dictated by the
615  // assembly map. We loop over each element, then the boundary degrees of
616  // freedom for u, boundary for v, followed by the interior for u and then
617  // interior for v.
618  // - inout holds the Dirichlet degrees of freedom in the local ordering,
619  // which have the boundary conditions imposed
620  Array<OneD, NekDouble> forCoeffs(nVel * nCoeffs, 0.0);
621  Array<OneD, NekDouble> inout (nVel * nCoeffs, 0.0);
622  Array<OneD, NekDouble> tmp (nCoeffs);
623 
624  for (nv = 0; nv < nVel; ++nv)
625  {
626  // Take the inner product of the forcing function.
627  m_fields[nv]->IProductWRTBase_IterPerExp(forcing[nv], tmp);
628 
629  // Impose Dirichlet condition on field which should be initialised
630  Array<OneD, NekDouble> loc_inout = m_fields[nv]->UpdateCoeffs();
631  m_fields[nv]->ImposeDirichletConditions(loc_inout);
632 
633  // Scatter forcing into RHS vector according to the ordering dictated in
634  // the comment above.
635  for (i = 0; i < m_fields[nv]->GetExpSize(); ++i)
636  {
637  int nBnd = m_bmap[i].size();
638  int nInt = m_imap[i].size();
639  int offset = m_fields[nv]->GetCoeff_Offset(i);
640 
641  for (j = 0; j < nBnd; ++j)
642  {
643  forCoeffs[nVel*offset + nv*nBnd + j] =
644  -1.0*tmp[offset+m_bmap[i][j]];
645  inout[nVel*offset + nv*nBnd + j] =
646  loc_inout[offset+m_bmap[i][j]];
647  }
648  for (j = 0; j < nInt; ++j)
649  {
650  forCoeffs[nVel*(offset + nBnd) + nv*nInt + j] =
651  -1.0*tmp[offset+m_imap[i][j]];
652  }
653  }
654  }
655 
656  //
657  // -- Perform solve
658  //
659  // Solve.
660  linSys->Solve(forCoeffs, inout, m_assemblyMap);
661 
662  //
663  // -- Postprocess
664  //
665 
666  // Scatter back to field degrees of freedom
667  for (nv = 0; nv < nVel; ++nv)
668  {
669  for (i = 0; i < m_fields[nv]->GetExpSize(); ++i)
670  {
671  int nBnd = m_bmap[i].size();
672  int nInt = m_imap[i].size();
673  int offset = m_fields[nv]->GetCoeff_Offset(i);
674 
675  for (j = 0; j < nBnd; ++j)
676  {
677  m_fields[nv]->UpdateCoeffs()[offset+m_bmap[i][j]] =
678  inout[nVel*offset + nv*nBnd + j];
679  }
680  for (j = 0; j < nInt; ++j)
681  {
682  m_fields[nv]->UpdateCoeffs()[offset+m_imap[i][j]] =
683  inout[nVel*(offset + nBnd) + nv*nInt + j];
684  }
685  }
686  m_fields[nv]->BwdTrans(m_fields[nv]->GetCoeffs(),
687  m_fields[nv]->UpdatePhys());
688  }
689 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
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:370
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
static PreconditionerSharedPtr NullPreconditionerSharedPtr
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:50
@ eDeformed
Geometry is curved or has non-constant factors.
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:51
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.cpp:322
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267

References ASSERTL0, 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 
)
protectedvirtual

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 823 of file LinearElasticSystem.cpp.

826 {
827  const int nVel = m_fields[0]->GetCoordim(0);
828  const int nCoeffs = m_fields[0]->GetNcoeffs();
829 
830  if (m_temperature.size() == 0)
831  {
832  return;
833  }
834 
835  for (int i = 0; i < nVel; ++i)
836  {
837  Array<OneD, NekDouble> tFwd(nCoeffs);
838  m_fields[i]->FwdTrans(m_temperature[i], tFwd);
839  fieldcoeffs.push_back(tFwd);
840  variables.push_back(
841  "ThermStressDiv" + boost::lexical_cast<std::string>(i));
842  }
843 
844  if (m_stress.size() == 0)
845  {
846  return;
847  }
848 
849  for (int i = 0; i < nVel; ++i)
850  {
851  for (int j = 0; j < nVel; ++j)
852  {
853  Array<OneD, NekDouble> tFwd(nCoeffs);
854  m_fields[i]->FwdTrans(m_stress[i][j], tFwd);
855  fieldcoeffs.push_back(tFwd);
856  variables.push_back(
857  "ThermStress"
858  + boost::lexical_cast<std::string>(i)
859  + boost::lexical_cast<std::string>(j));
860  }
861  }
862 }

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

◆ v_GenerateSummary()

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

Generate summary at runtime.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Reimplemented in Nektar::IterativeElasticSystem.

Definition at line 352 of file LinearElasticSystem.cpp.

353 {
355 
356  AddSummaryItem(s, "Young's modulus", m_E);
357  AddSummaryItem(s, "Poisson ratio", m_nu);
358 }
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 ( )
protectedvirtual

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.

Reimplemented in Nektar::IterativeElasticSystem.

Definition at line 140 of file LinearElasticSystem.cpp.

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

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

string Nektar::LinearElasticSystem::className
static
Initial value:
RegisterCreatorFunction("LinearElasticSystem",
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 67 of file LinearElasticSystem.h.

◆ m_assemblyMap

CoupledAssemblyMapSharedPtr Nektar::LinearElasticSystem::m_assemblyMap
protected

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

Definition at line 94 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 92 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 98 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 104 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 100 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 102 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 90 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 106 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 88 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 96 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 110 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 108 of file LinearElasticSystem.h.

Referenced by v_DoSolve(), and v_ExtraFldOutput().