Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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:
Inheritance graph
[legend]
Collaboration diagram for Nektar::LinearElasticSystem:
Collaboration graph
[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 >
boost::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 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 void EvaluateFunction (Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
 Evaluates a function as specified in the session file. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::vector< std::string > pFieldNames, Array< OneD, Array< OneD, NekDouble > > &pFields, const std::string &pName, const NekDouble &pTime=0.0, const int domain=0)
 Populate given fields with the function from session. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::vector< std::string > pFieldNames, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const std::string &pName, const NekDouble &pTime=0.0, const int domain=0)
 Populate given fields with the function from session. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, const NekDouble &pTime=0.0, const int domain=0)
 
SOLVER_UTILS_EXPORT std::string DescribeFunction (std::string pFieldName, const std::string &pFunctionName, const int domain)
 Provide a description of a function for a given field name. More...
 
SOLVER_UTILS_EXPORT void InitialiseBaseFlow (Array< OneD, Array< OneD, NekDouble > > &base)
 Perform initialisation of the base flow. 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, NekDouble
ErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionGreensDivergenceForm (const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
 Compute the inner product $ (\nabla \phi \cdot F) $. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionDivergenceForm (const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
 Compute the inner product $ (\phi, \nabla \cdot F) $. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionNonConservativeForm (const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, bool UseContCoeffs=false)
 Compute the inner product $ (\phi, V\cdot \nabla u) $. More...
 
f SOLVER_UTILS_EXPORT void AdvectionNonConservativeForm (const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
 Compute the non-conservative advection. More...
 
SOLVER_UTILS_EXPORT void WeakDGAdvection (const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false, int nvariables=0)
 Calculate the weak discontinuous Galerkin advection. More...
 
SOLVER_UTILS_EXPORT void WeakDGDiffusion (const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false)
 Calculate weak DG Diffusion in the LDG form. 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 ScanForHistoryPoints ()
 Builds map of which element holds each history point. More...
 
SOLVER_UTILS_EXPORT void WriteHistoryData (std::ostream &out)
 Probe each history point and write to 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::FieldMetaDataMap
UpdateFieldMetaDataMap ()
 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 GetNumElmVelocity ()
 
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 GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &fluxX, Array< OneD, Array< OneD, NekDouble > > &fluxY)
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, const int j, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
 
SOLVER_UTILS_EXPORT void NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
 
SOLVER_UTILS_EXPORT void NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxX, Array< OneD, Array< OneD, NekDouble > > &numfluxY)
 
SOLVER_UTILS_EXPORT void NumFluxforScalar (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
 
SOLVER_UTILS_EXPORT void NumFluxforVector (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int NoCaseStringCompare (const std::string &s1, const std::string &s2)
 Perform a case-insensitive string comparison. More...
 
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)
 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)
 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)
 Initialises EquationSystem class members. More...
 
int nocase_cmp (const std::string &s1, const std::string &s2)
 
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)
 
SOLVER_UTILS_EXPORT void SetUpBaseFields (SpatialDomains::MeshGraphSharedPtr &mesh)
 
SOLVER_UTILS_EXPORT void ImportFldBase (std::string pInfile, SpatialDomains::MeshGraphSharedPtr pGraph)
 
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...
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
std::map< std::string, Array
< OneD, Array< OneD, float > > > 
m_interpWeights
 Map of the interpolation weights for a specific filename. More...
 
std::map< std::string, Array
< OneD, Array< OneD, unsigned
int > > > 
m_interpInds
 Map of the interpolation indices for a specific filename. More...
 
Array< OneD,
MultiRegions::ExpListSharedPtr
m_fields
 Array holding all dependent variables. More...
 
Array< OneD,
MultiRegions::ExpListSharedPtr
m_base
 Base fields. More...
 
Array< OneD,
MultiRegions::ExpListSharedPtr
m_derivedfields
 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...
 
std::set< std::string > m_loadedFields
 
NekDouble m_checktime
 Time between checkpoints. More...
 
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, Array< OneD,
Array< OneD, NekDouble > > > 
m_gradtan
 1 x nvariable x nq More...
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 
m_tanbasis
 2 x m_spacedim x nq 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...
 

Detailed Description

Base class for linear elastic system.

Definition at line 50 of file LinearElasticSystem.h.

Constructor & Destructor Documentation

Nektar::LinearElasticSystem::LinearElasticSystem ( const LibUtilities::SessionReaderSharedPtr pSession)
protected

Default constructor.

Definition at line 129 of file LinearElasticSystem.cpp.

131  : EquationSystem(pSession)
132 {
133 }
SOLVER_UTILS_EXPORT EquationSystem(const LibUtilities::SessionReaderSharedPtr &pSession)
Initialises EquationSystem class members.

Member Function Documentation

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

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

Referenced by BuildMatrixSystem().

885 {
886  const int nCoeffs = exp->GetNcoeffs();
887  const int nPhys = exp->GetTotPoints();
888  int i;
889 
891  nCoeffs, nCoeffs, 0.0, eFULL);
892 
893  Array<OneD, NekDouble> tmp2(nPhys);
894  Array<OneD, NekDouble> tmp3(nPhys);
895 
896  for (i = 0; i < nCoeffs; ++i)
897  {
898  Array<OneD, NekDouble> tmp1(nCoeffs, 0.0);
899  tmp1[i] = 1.0;
900 
901  exp->BwdTrans ( tmp1, tmp2);
902  exp->PhysDeriv (k1, tmp2, tmp3);
903  exp->IProductWRTDerivBase(k2, tmp3, tmp1);
904 
905  Vmath::Smul(
906  nCoeffs, scale, &tmp1[0], 1, &(ret->GetPtr())[0]+i*nCoeffs, 1);
907  }
908 
909  return ret;
910 }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
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 240 of file LinearElasticSystem.cpp.

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().

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

Creates an instance of this class.

Definition at line 57 of file LinearElasticSystem.h.

59  {
60  EquationSystemSharedPtr p = MemoryManager<
61  LinearElasticSystem>::AllocateSharedPtr(pSession);
62  p->InitObject();
63  return p;
64  }
boost::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
LinearElasticSystem(const LibUtilities::SessionReaderSharedPtr &pSession)
Default constructor.
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 792 of file LinearElasticSystem.cpp.

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

Referenced by BuildMatrixSystem().

796 {
797  int i, j, k, l;
798  const int nVel = mat.GetRows();
799  const int nB = exp->NumBndryCoeffs();
800  const int nI = exp->GetNcoeffs() - nB;
801  const int nBnd = exp->NumBndryCoeffs() * nVel;
802  const int nInt = exp->GetNcoeffs() * nVel - nBnd;
803  const MatrixStorage s = eFULL; // Maybe look into doing symmetric
804  // version of this?
805 
806  DNekMatSharedPtr A =
807  MemoryManager<DNekMat>::AllocateSharedPtr(nBnd, nBnd, 0.0, s);
808  DNekMatSharedPtr B =
809  MemoryManager<DNekMat>::AllocateSharedPtr(nBnd, nInt, 0.0, s);
810  DNekMatSharedPtr C =
811  MemoryManager<DNekMat>::AllocateSharedPtr(nInt, nBnd, 0.0, s);
812  DNekMatSharedPtr D =
813  MemoryManager<DNekMat>::AllocateSharedPtr(nInt, nInt, 0.0, s);
814 
815  for (i = 0; i < nVel; ++i)
816  {
817  for (j = 0; j < nVel; ++j)
818  {
819  // Boundary-boundary and boundary-interior
820  for (k = 0; k < nB; ++k)
821  {
822  for (l = 0; l < nB; ++l)
823  {
824  (*A)(k + i*nB, l + j*nB) =
825  (*mat[i][j])(m_bmap[n][k], m_bmap[n][l]);
826  }
827 
828  for (l = 0; l < nI; ++l)
829  {
830  (*B)(k + i*nB, l + j*nI) =
831  (*mat[i][j])(m_bmap[n][k], m_imap[n][l]);
832  }
833  }
834 
835  // Interior-boundary / interior-interior
836  for (k = 0; k < nI; ++k)
837  {
838  for (l = 0; l < nB; ++l)
839  {
840  (*C)(k + i*nI, l + j*nB) =
841  (*mat[i][j])(m_imap[n][k], m_bmap[n][l]);
842  }
843 
844  for (l = 0; l < nI; ++l)
845  {
846  (*D)(k + i*nI, l + j*nI) =
847  (*mat[i][j])(m_imap[n][k], m_imap[n][l]);
848  }
849  }
850  }
851  }
852 
853  // Construct static condensation matrices.
854  D->Invert();
855  (*B) = (*B)*(*D);
856  (*A) = (*A) - (*B)*(*C);
857 
858  DNekScalMatSharedPtr tmp_mat;
859  m_schurCompl->SetBlock(
860  n, n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(
861  1.0, A));
862  m_BinvD ->SetBlock(
863  n, n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(
864  1.0, B));
865  m_C ->SetBlock(
866  n, n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(
867  1.0, C));
868  m_Dinv ->SetBlock(
869  n, n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(
870  1.0, D));
871 }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
DNekScalBlkMatSharedPtr m_BinvD
Matrix of elemental components.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
DNekScalBlkMatSharedPtr m_Dinv
Matrix of elemental components.
Array< OneD, Array< OneD, unsigned int > > m_imap
Interior maps for each of the fields.
DNekScalBlkMatSharedPtr m_schurCompl
Schur complement boundary-boundary matrix.
DNekScalBlkMatSharedPtr m_C
Matrix of elemental components.
Array< OneD, Array< OneD, unsigned int > > m_bmap
Boundary maps for each of the fields.
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 391 of file LinearElasticSystem.cpp.

References ASSERTL0, BuildMatrixSystem(), Nektar::SpatialDomains::eDeformed, Nektar::eDIAGONAL, Nektar::MultiRegions::eDirectStaticCond, Nektar::SpatialDomains::eDirichlet, Nektar::eFULL, Nektar::MultiRegions::eIterativeStaticCond, Nektar::StdRegions::eLinearAdvectionReaction, Nektar::MultiRegions::ePETScStaticCond, Nektar::SolverUtils::EquationSystem::EvaluateFunction(), Nektar::MultiRegions::eXxtStaticCond, Vmath::Fill(), Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::iterator, 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(), Vmath::Neg(), Nektar::MultiRegions::NullPreconditionerSharedPtr, sign, Vmath::Smul(), Nektar::Transpose(), Vmath::Vadd(), and Vmath::Vcopy().

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

392 {
393 
394  int i, j, k, l, nv;
395  const int nVel = m_fields[0]->GetCoordim(0);
396 
397  // Build initial matrix system.
399 
400  // Now we've got the matrix system set up, create a GlobalLinSys
401  // object. We mask ourselves as LinearAdvectionReaction to create a full
402  // matrix instead of symmetric storage.
403  MultiRegions::GlobalLinSysKey key(
406 
407  // Currently either direct or iterative static condensation is
408  // supported.
409  if (m_assemblyMap->GetGlobalSysSolnType() == MultiRegions::eDirectStaticCond)
410  {
411  linSys = MemoryManager<
412  MultiRegions::GlobalLinSysDirectStaticCond>::AllocateSharedPtr(
413  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
414  m_assemblyMap);
415  }
416  else if (m_assemblyMap->GetGlobalSysSolnType() ==
418  {
419  linSys = MemoryManager<
420  MultiRegions::GlobalLinSysIterativeStaticCond>::AllocateSharedPtr(
421  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
423  }
424 #ifdef NEKTAR_USE_PETSC
425  else if (m_assemblyMap->GetGlobalSysSolnType() ==
427  {
428  linSys = MemoryManager<
429  MultiRegions::GlobalLinSysPETScStaticCond>::AllocateSharedPtr(
430  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
431  m_assemblyMap);
432  }
433 #endif
434 #ifdef NEKTAR_USE_MPI
435  else if (m_assemblyMap->GetGlobalSysSolnType() ==
437  {
438  linSys = MemoryManager<
439  MultiRegions::GlobalLinSysXxtStaticCond>::AllocateSharedPtr(
440  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
441  m_assemblyMap);
442  }
443 #endif
444 
445  linSys->Initialise(m_assemblyMap);
446 
447  const int nCoeffs = m_fields[0]->GetNcoeffs();
448  const int nGlobDofs = m_assemblyMap->GetNumGlobalCoeffs();
449 
450  //
451  // -- Evaluate forcing functions
452  //
453 
454  // Evaluate the forcing function from the XML file.
455  Array<OneD, Array<OneD, NekDouble> > forcing(nVel);
456  EvaluateFunction(forcing, "Forcing");
457 
458  // Add temperature term
459  string tempEval;
460  m_session->LoadSolverInfo("Temperature", tempEval, "None");
461 
462  if (tempEval == "Jacobian")
463  {
464  // Allocate storage
465  m_temperature = Array<OneD, Array<OneD, NekDouble> >(nVel);
466 
467  for (nv = 0; nv < nVel; ++nv)
468  {
469  Array<OneD, NekDouble> tmp;
470  m_temperature[nv] = Array<OneD, NekDouble>(
471  m_fields[nv]->GetNpoints());
472 
473  for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
474  {
475  // Calculate element area
477  m_fields[0]->GetExp(i);
479  exp->GetPointsKeys();
480  Array<OneD, NekDouble> jac =
481  exp->GetMetricInfo()->GetJac(pkey);
482 
483  int offset = m_fields[0]->GetPhys_Offset(i);
484 
485  if (exp->GetMetricInfo()->GetGtype() ==
487  {
488  Vmath::Smul(exp->GetTotPoints(), m_beta, jac, 1,
489  tmp = m_temperature[nv] + offset, 1);
490  }
491  else
492  {
493  Vmath::Fill(exp->GetTotPoints(), m_beta*jac[0],
494  tmp = m_temperature[nv] + offset, 1);
495  }
496  }
497  m_fields[nv]->PhysDeriv(nv, m_temperature[nv], forcing[nv]);
498  }
499  }
500  else if (tempEval == "Metric")
501  {
502  ASSERTL0((m_fields[0]->GetCoordim(0) == 2 &&
503  m_graph->GetAllQuadGeoms().size() == 0) ||
504  (m_fields[0]->GetCoordim(0) == 3 &&
505  m_graph->GetAllPrismGeoms().size() == 0 &&
506  m_graph->GetAllPyrGeoms ().size() == 0 &&
507  m_graph->GetAllHexGeoms ().size() == 0),
508  "LinearIdealMetric temperature only implemented for "
509  "two-dimensional triangular meshes or three-dimensional "
510  "tetrahedral meshes.");
511 
512  m_temperature = Array<OneD, Array<OneD, NekDouble> >(nVel);
513  m_stress = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >(nVel);
514 
515  for (nv = 0; nv < nVel; ++nv)
516  {
517  m_temperature[nv] = Array<OneD, NekDouble>(
518  m_fields[nv]->GetNpoints(), 0.0);
519  m_stress[nv] = Array<OneD, Array<OneD, NekDouble> >(nVel);
520 
521  for (i = 0; i < nVel; ++i)
522  {
523  m_stress[nv][i] = Array<OneD, NekDouble>(
524  m_fields[nv]->GetNpoints(), 0.0);
525  }
526  }
527 
528  for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
529  {
530  // Calculate element area
532  m_fields[0]->GetExp(i);
533  LibUtilities::PointsKeyVector pkey = exp->GetPointsKeys();
534  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > deriv =
535  exp->GetMetricInfo()->GetDeriv(pkey);
536  int offset = m_fields[0]->GetPhys_Offset(i);
537 
538  DNekMat i2rm = MappingIdealToRef(exp->GetGeom());
539 
540  // Compute metric tensor
541  DNekMat jac (nVel, nVel, 0.0, eFULL);
542  DNekMat jacIdeal(nVel, nVel, 0.0, eFULL);
543  DNekMat metric (nVel, nVel, 0.0, eFULL);
544 
545  for (j = 0; j < deriv[0][0].num_elements(); ++j)
546  {
547  for (k = 0; k < nVel; ++k)
548  {
549  for (l = 0; l < nVel; ++l)
550  {
551  jac(l,k) = deriv[k][l][j];
552  }
553  }
554 
555  jacIdeal = jac * i2rm;
556  metric = Transpose(jacIdeal) * jacIdeal;
557 
558  // Compute eigenvalues/eigenvectors of metric tensor using
559  // ideal mapping.
560  char jobvl = 'N', jobvr = 'V';
561  int worklen = 8*nVel, info;
562 
563  DNekMat eval (nVel, nVel, 0.0, eDIAGONAL);
564  DNekMat evec (nVel, nVel, 0.0, eFULL);
565  DNekMat evecinv(nVel, nVel, 0.0, eFULL);
566  Array<OneD, NekDouble> vl (nVel*nVel);
567  Array<OneD, NekDouble> work(worklen);
568  Array<OneD, NekDouble> wi (nVel);
569 
570  Lapack::Dgeev(jobvl, jobvr, nVel, metric.GetRawPtr(), nVel,
571  &(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
572  &(evec.GetPtr())[0], nVel,
573  &work[0], worklen, info);
574 
575  evecinv = evec;
576  evecinv.Invert();
577 
578  // rescaling of the eigenvalues
579  for (nv = 0; nv < nVel; ++nv)
580  {
581  eval(nv,nv) = m_beta * (sqrt(eval(nv,nv)) - 1.0);
582  }
583 
584  DNekMat beta = evec * eval * evecinv;
585 
586  for (k = 0; k < nVel; ++k)
587  {
588  for (l = 0; l < nVel; ++l)
589  {
590  m_stress[k][l][offset+j] = beta(k,l);
591  }
592  }
593  }
594 
595  if (deriv[0][0].num_elements() != exp->GetTotPoints())
596  {
597  Array<OneD, NekDouble> tmp;
598  for (k = 0; k < nVel; ++k)
599  {
600  for (l = 0; l < nVel; ++l)
601  {
602  Vmath::Fill(
603  exp->GetTotPoints(), m_stress[k][l][offset],
604  tmp = m_stress[k][l] + offset, 1);
605  }
606  }
607  }
608  }
609 
610  // Calculate divergence of stress tensor.
611  Array<OneD, NekDouble> tmpderiv(m_fields[0]->GetNpoints());
612  for (i = 0; i < nVel; ++i)
613  {
614  for (j = 0; j < nVel; ++j)
615  {
616  m_fields[i]->PhysDeriv(j, m_stress[i][j], tmpderiv);
617  Vmath::Vadd (m_fields[i]->GetNpoints(), tmpderiv, 1,
618  m_temperature[i], 1, m_temperature[i], 1);
619  }
620 
622  m_temperature[i], 1, forcing[i], 1);
623  }
624  }
625  else if (tempEval != "None")
626  {
627  ASSERTL0(false, "Unknown temperature form: " + tempEval);
628  }
629 
630  // Set up some temporary storage.
631  //
632  // - forCoeffs holds the forcing coefficients in a local ordering;
633  // however note that the ordering is different and dictated by the
634  // assembly map. We loop over each element, then the boundary degrees of
635  // freedom for u, boundary for v, followed by the interior for u and then
636  // interior for v.
637  // - rhs is the global assembly of forCoeffs.
638  // - inout holds the Dirichlet degrees of freedom in the global ordering,
639  // which have been assembled from the boundary expansion.
640  Array<OneD, NekDouble> forCoeffs(nVel * nCoeffs, 0.0);
641  Array<OneD, NekDouble> inout (nGlobDofs, 0.0);
642  Array<OneD, NekDouble> rhs (nGlobDofs, 0.0);
643 
644  for (nv = 0; nv < nVel; ++nv)
645  {
646  // Take the inner product of the forcing function.
647  Array<OneD, NekDouble> tmp(nCoeffs);
648  m_fields[nv]->IProductWRTBase_IterPerExp(forcing[nv], tmp);
649 
650  // Scatter forcing into RHS vector according to the ordering dictated in
651  // the comment above.
652  for (i = 0; i < m_fields[nv]->GetExpSize(); ++i)
653  {
654  int nBnd = m_bmap[i].num_elements();
655  int nInt = m_imap[i].num_elements();
656  int offset = m_fields[nv]->GetCoeff_Offset(i);
657 
658  for (j = 0; j < nBnd; ++j)
659  {
660  forCoeffs[nVel*offset + nv*nBnd + j] =
661  tmp[offset+m_bmap[i][j]];
662  }
663  for (j = 0; j < nInt; ++j)
664  {
665  forCoeffs[nVel*(offset + nBnd) + nv*nInt + j] =
666  tmp[offset+m_imap[i][j]];
667  }
668  }
669  }
670 
671  // -- Impose Dirichlet boundary conditions.
672 
673  // First try to do parallel assembly: the intention here is that Dirichlet
674  // values at some edges/vertices need to be communicated to processes which
675  // don't contain the entire boundary region. See
676  // ContField2D::v_ImposeDirichletConditions for more detail.
677  map<int, vector<MultiRegions::ExtraDirDof> > &extraDirDofs =
678  m_assemblyMap->GetExtraDirDofs();
679  map<int, vector<MultiRegions::ExtraDirDof> >::iterator it;
680 
681  for (nv = 0; nv < nVel; ++nv)
682  {
683  const Array<OneD, const MultiRegions::ExpListSharedPtr> &bndCondExp
684  = m_fields[nv]->GetBndCondExpansions();
685 
686  // First try to do parallel stuff
687  for (it = extraDirDofs.begin(); it != extraDirDofs.end(); ++it)
688  {
689  for (i = 0; i < it->second.size(); ++i)
690  {
691  inout[it->second.at(i).get<1>()*nVel + nv] =
692  bndCondExp[it->first]->GetCoeffs()[
693  it->second.at(i).get<0>()]*it->second.at(i).get<2>();
694  }
695  }
696  }
697 
698  m_assemblyMap->UniversalAssemble(inout);
699 
700  // Counter for the local Dirichlet boundary to global ordering.
701  int bndcnt = 0;
702 
703  // Now assemble local boundary contributions.
704  for (nv = 0; nv < nVel; ++nv)
705  {
706  const Array<OneD, const MultiRegions::ExpListSharedPtr> &bndCondExp
707  = m_fields[nv]->GetBndCondExpansions();
708  const Array<OneD, const int> &bndMap
709  = m_assemblyMap->GetBndCondCoeffsToGlobalCoeffsMap();
710 
711  for (i = 0; i < bndCondExp.num_elements(); ++i)
712  {
713  if (m_fields[nv]->GetBndConditions()[i]->GetBoundaryConditionType()
715  {
716  const Array<OneD,const NekDouble> &bndCoeffs =
717  bndCondExp[i]->GetCoeffs();
718 
719  for (j = 0; j < bndCondExp[i]->GetNcoeffs(); ++j)
720  {
721  NekDouble sign =
722  m_assemblyMap->GetBndCondCoeffsToGlobalCoeffsSign(
723  bndcnt);
724  inout[bndMap[bndcnt++]] = sign * bndCoeffs[j];
725  }
726  }
727  else
728  {
729  bndcnt += bndCondExp[i]->GetNcoeffs();
730  }
731  }
732  }
733 
734  //
735  // -- Perform solve
736  //
737 
738  // Assemble forcing into the RHS.
739  m_assemblyMap->Assemble(forCoeffs, rhs);
740 
741  // Negate RHS to be consistent with matrix definition.
742  Vmath::Neg(rhs.num_elements(), rhs, 1);
743 
744  // Solve.
745  linSys->Solve(rhs, inout, m_assemblyMap);
746 
747  // Scatter the global ordering back to the alternate local ordering.
748  Array<OneD, NekDouble> tmp(nVel * nCoeffs);
749  m_assemblyMap->GlobalToLocal(inout, tmp);
750 
751  //
752  // -- Postprocess
753  //
754 
755  // Scatter back to field degrees of freedom
756  for (nv = 0; nv < nVel; ++nv)
757  {
758  for (i = 0; i < m_fields[nv]->GetExpSize(); ++i)
759  {
760  int nBnd = m_bmap[i].num_elements();
761  int nInt = m_imap[i].num_elements();
762  int offset = m_fields[nv]->GetCoeff_Offset(i);
763 
764  for (j = 0; j < nBnd; ++j)
765  {
766  m_fields[nv]->UpdateCoeffs()[offset+m_bmap[i][j]] =
767  tmp[nVel*offset + nv*nBnd + j];
768  }
769  for (j = 0; j < nInt; ++j)
770  {
771  m_fields[nv]->UpdateCoeffs()[offset+m_imap[i][j]] =
772  tmp[nVel*(offset + nBnd) + nv*nInt + j];
773  }
774  }
775  m_fields[nv]->BwdTrans(m_fields[nv]->GetCoeffs(),
776  m_fields[nv]->UpdatePhys());
777  }
778 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
NekDouble m_beta
Parameter dictating amount of thermal stress to add.
DNekScalBlkMatSharedPtr m_BinvD
Matrix of elemental components.
Array< OneD, Array< OneD, NekDouble > > m_temperature
Storage for the temperature terms.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_stress
Storage for the thermal stress terms.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
static PreconditionerSharedPtr NullPreconditionerSharedPtr
DNekScalBlkMatSharedPtr m_Dinv
Matrix of elemental components.
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:52
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Array< OneD, Array< OneD, unsigned int > > m_imap
Interior maps for each of the fields.
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
DNekScalBlkMatSharedPtr m_schurCompl
Schur complement boundary-boundary matrix.
DNekScalBlkMatSharedPtr m_C
Matrix of elemental components.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
CoupledAssemblyMapSharedPtr m_assemblyMap
Assembly map for the coupled (u,v,w) system.
Array< OneD, Array< OneD, unsigned int > > m_bmap
Boundary maps for each of the fields.
boost::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:52
void BuildMatrixSystem()
Build matrix system for linear elasticity equations.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
DNekMat MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Geometry is curved or has non-constant factors.
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:285
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 912 of file LinearElasticSystem.cpp.

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

915 {
916  const int nVel = m_fields[0]->GetCoordim(0);
917  const int nCoeffs = m_fields[0]->GetNcoeffs();
918 
919  if (m_temperature.num_elements() == 0)
920  {
921  return;
922  }
923 
924  for (int i = 0; i < nVel; ++i)
925  {
926  Array<OneD, NekDouble> tFwd(nCoeffs);
927  m_fields[i]->FwdTrans(m_temperature[i], tFwd);
928  fieldcoeffs.push_back(tFwd);
929  variables.push_back(
930  "ThermStressDiv" + boost::lexical_cast<std::string>(i));
931  }
932 
933  if (m_stress.num_elements() == 0)
934  {
935  return;
936  }
937 
938  for (int i = 0; i < nVel; ++i)
939  {
940  for (int j = 0; j < nVel; ++j)
941  {
942  Array<OneD, NekDouble> tFwd(nCoeffs);
943  m_fields[i]->FwdTrans(m_stress[i][j], tFwd);
944  fieldcoeffs.push_back(tFwd);
945  variables.push_back(
946  "ThermStress"
947  + boost::lexical_cast<std::string>(i)
948  + boost::lexical_cast<std::string>(j));
949  }
950  }
951 }
Array< OneD, Array< OneD, NekDouble > > m_temperature
Storage for the temperature terms.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_stress
Storage for the thermal stress terms.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
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 369 of file LinearElasticSystem.cpp.

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

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

370 {
372 
373  AddSummaryItem(s, "Young's modulus", m_E);
374  AddSummaryItem(s, "Poisson ratio", m_nu);
375 }
NekDouble m_nu
Poisson ratio.
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:50
NekDouble m_E
Young's modulus.
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
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 142 of file LinearElasticSystem.cpp.

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().

143 {
145 
146  const int nVel = m_fields[0]->GetCoordim(0);
147  int n;
148 
149  ASSERTL0(nVel > 1, "Linear elastic solver not set up for"
150  " this dimension (only 2D/3D supported).");
151 
152  // Make sure that we have Young's modulus and Poisson ratio set.
153  m_session->LoadParameter("E", m_E, 1.00);
154  m_session->LoadParameter("nu", m_nu, 0.25);
155  m_session->LoadParameter("Beta", m_beta, 1.00);
156 
157  // Create a coupled assembly map which allows us to tie u, v and w
158  // fields together.
159  if (nVel == 2)
160  {
161  MultiRegions::ContField2DSharedPtr u = boost::dynamic_pointer_cast<
162  MultiRegions::ContField2D>(m_fields[0]);
165  m_graph,
166  u->GetLocalToGlobalMap(),
167  m_fields[0]->GetBndConditions(),
168  m_fields);
169  }
170 
171  if (nVel == 3)
172  {
173  MultiRegions::ContField3DSharedPtr u = boost::dynamic_pointer_cast<
174  MultiRegions::ContField3D>(m_fields[0]);
175  m_assemblyMap = MemoryManager<CoupledAssemblyMap>
177  m_graph,
178  u->GetLocalToGlobalMap(),
179  m_fields[0]->GetBndConditions(),
180  m_fields);
181  }
182 
183  // Figure out size of our new matrix systems by looping over all expansions
184  // and multiply number of coefficients by velocity components.
185  const int nEl = m_fields[0]->GetExpSize();
187 
188  Array<OneD, unsigned int> sizeBnd(nEl);
189  Array<OneD, unsigned int> sizeInt(nEl);
190 
191  // Allocate storage for boundary and interior map caches.
192  m_bmap = Array<OneD, Array<OneD, unsigned int> >(nEl);
193  m_imap = Array<OneD, Array<OneD, unsigned int> >(nEl);
194 
195  // Cache interior and boundary maps to avoid recalculating
196  // constantly. Should really be handled by manager in StdRegions but not
197  // really working yet.
198  for (n = 0; n < nEl; ++n)
199  {
200  exp = m_fields[0]->GetExp(m_fields[0]->GetOffset_Elmt_Id(n));
201  sizeBnd[n] = nVel * exp->NumBndryCoeffs();
202  sizeInt[n] = nVel * exp->GetNcoeffs() - sizeBnd[n];
203  exp->GetBoundaryMap(m_bmap[n]);
204  exp->GetInteriorMap(m_imap[n]);
205  }
206 
207  // Create block matrix storage for the statically condensed system.
208  MatrixStorage blkmatStorage = eDIAGONAL;
210  sizeBnd, sizeBnd, blkmatStorage);
212  sizeBnd, sizeInt, blkmatStorage);
214  sizeInt, sizeBnd, blkmatStorage);
216  sizeInt, sizeInt, blkmatStorage);
217 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
NekDouble m_nu
Poisson ratio.
NekDouble m_beta
Parameter dictating amount of thermal stress to add.
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:293
DNekScalBlkMatSharedPtr m_BinvD
Matrix of elemental components.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Initialisation object for EquationSystem.
DNekScalBlkMatSharedPtr m_Dinv
Matrix of elemental components.
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Array< OneD, Array< OneD, unsigned int > > m_imap
Interior maps for each of the fields.
DNekScalBlkMatSharedPtr m_schurCompl
Schur complement boundary-boundary matrix.
DNekScalBlkMatSharedPtr m_C
Matrix of elemental components.
NekDouble m_E
Young's modulus.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
CoupledAssemblyMapSharedPtr m_assemblyMap
Assembly map for the coupled (u,v,w) system.
Array< OneD, Array< OneD, unsigned int > > m_bmap
Boundary maps for each of the fields.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
boost::shared_ptr< ContField3D > ContField3DSharedPtr
Definition: ContField3D.h:191

Friends And Related Function Documentation

friend class MemoryManager< LinearElasticSystem >
friend

Class may only be instantiated through the MemoryManager.

Definition at line 54 of file LinearElasticSystem.h.

Member Data Documentation

string Nektar::LinearElasticSystem::className
static
Initial value:
RegisterCreatorFunction("LinearElasticSystem",

Name of class.

Definition at line 67 of file LinearElasticSystem.h.

CoupledAssemblyMapSharedPtr Nektar::LinearElasticSystem::m_assemblyMap
protected

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

Definition at line 93 of file LinearElasticSystem.h.

Referenced by v_DoSolve(), and v_InitObject().

NekDouble Nektar::LinearElasticSystem::m_beta
protected

Parameter dictating amount of thermal stress to add.

Definition at line 91 of file LinearElasticSystem.h.

Referenced by v_DoSolve(), and v_InitObject().

DNekScalBlkMatSharedPtr Nektar::LinearElasticSystem::m_BinvD
protected

Matrix of elemental $ B^{-1}D $ components.

Definition at line 97 of file LinearElasticSystem.h.

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

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

Boundary maps for each of the fields.

Definition at line 103 of file LinearElasticSystem.h.

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

DNekScalBlkMatSharedPtr Nektar::LinearElasticSystem::m_C
protected

Matrix of elemental $ C $ components.

Definition at line 99 of file LinearElasticSystem.h.

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

DNekScalBlkMatSharedPtr Nektar::LinearElasticSystem::m_Dinv
protected

Matrix of elemental $ D^{-1} $ components.

Definition at line 101 of file LinearElasticSystem.h.

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

NekDouble Nektar::LinearElasticSystem::m_E
protected

Young's modulus.

Definition at line 89 of file LinearElasticSystem.h.

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

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

Interior maps for each of the fields.

Definition at line 105 of file LinearElasticSystem.h.

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

NekDouble Nektar::LinearElasticSystem::m_nu
protected

Poisson ratio.

Definition at line 87 of file LinearElasticSystem.h.

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

DNekScalBlkMatSharedPtr Nektar::LinearElasticSystem::m_schurCompl
protected

Schur complement boundary-boundary matrix.

Definition at line 95 of file LinearElasticSystem.h.

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

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

Storage for the thermal stress terms.

Definition at line 109 of file LinearElasticSystem.h.

Referenced by v_DoSolve(), and v_ExtraFldOutput().

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

Storage for the temperature terms.

Definition at line 107 of file LinearElasticSystem.h.

Referenced by v_DoSolve(), and v_ExtraFldOutput().