55 #ifdef NEKTAR_USE_PETSC
65 RegisterCreatorFunction(
"LinearElasticSystem",
66 LinearElasticSystem::create);
80 int n = geom->GetNumVerts(), i, j;
86 for (i = 0; i < n; ++i)
88 for (j = 0; j < n-1; ++j)
90 map(j,i) = (*geom->GetVertex(i))[j];
97 mapref(0,0) = -1.0; mapref(1,0) = -1.0;
98 mapref(0,1) = 1.0; mapref(1,1) = -1.0;
99 mapref(0,2) = -1.0; mapref(1,2) = 1.0;
103 mapref(0,0) = -1.0; mapref(1,0) = -1.0; mapref(2,0) = -1.0;
104 mapref(0,1) = 1.0; mapref(1,1) = -1.0; mapref(2,1) = -1.0;
105 mapref(0,2) = -1.0; mapref(1,2) = 1.0; mapref(2,2) = -1.0;
106 mapref(0,3) = -1.0; mapref(1,3) = -1.0; mapref(2,3) = 1.0;
114 for (i = 0; i < n-1; ++i)
116 for (j = 0; j < n-1; ++j)
118 mapred(i,j) = newmap(i,j);
129 LinearElasticSystem::LinearElasticSystem(
146 const int nVel =
m_fields[0]->GetCoordim(0);
149 ASSERTL0(nVel > 1,
"Linear elastic solver not set up for"
150 " this dimension (only 2D/3D supported).");
166 u->GetLocalToGlobalMap(),
178 u->GetLocalToGlobalMap(),
185 const int nEl =
m_fields[0]->GetExpSize();
198 for (n = 0; n < nEl; ++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]);
210 sizeBnd, sizeBnd, blkmatStorage);
212 sizeBnd, sizeInt, blkmatStorage);
214 sizeInt, sizeBnd, blkmatStorage);
216 sizeInt, sizeInt, blkmatStorage);
242 const int nEl =
m_fields[0]->GetExpSize();
243 const int nVel =
m_fields[0]->GetCoordim(0);
255 bool verbose =
m_session->DefinesCmdLineArgument(
"verbose");
256 bool root =
m_session->GetComm()->GetRank() == 0;
261 for (n = 0; n < nEl; ++n)
264 const int nPhys = exp->GetTotPoints();
276 *exp, factors, varcoeffA);
279 *exp, factors, varcoeffD);
286 mat[0][0] = exp->GenMatrix(matkeyA);
288 mat[1][0] = mat[0][1];
289 mat[1][1] = exp->GenMatrix(matkeyD);
296 cout <<
"\rBuilding matrix system: "
297 << (int)(100.0 * n / nEl) <<
"%" << flush;
303 for (n = 0; n < nEl; ++n)
306 const int nPhys = exp->GetTotPoints();
323 *exp, factors, varcoeffA);
326 *exp, factors, varcoeffE);
329 *exp, factors, varcoeffI);
337 mat[0][0] = exp->GenMatrix(matkeyA);
341 mat[1][0] = mat[0][1];
342 mat[1][1] = exp->GenMatrix(matkeyE);
345 mat[2][0] = mat[0][2];
346 mat[2][1] = mat[1][2];
347 mat[2][2] = exp->GenMatrix(matkeyI);
354 cout <<
"\rBuilding matrix system: "
355 << (int)(100.0 * n / nEl) <<
"%" << flush;
362 cout <<
"\rBuilding matrix system: done." << endl;
395 const int nVel =
m_fields[0]->GetCoordim(0);
424 #ifdef NEKTAR_USE_PETSC
434 #ifdef NEKTAR_USE_MPI
447 const int nCoeffs =
m_fields[0]->GetNcoeffs();
460 m_session->LoadSolverInfo(
"Temperature", tempEval,
"None");
462 if (tempEval ==
"Jacobian")
467 for (nv = 0; nv < nVel; ++nv)
473 for (i = 0; i <
m_fields[0]->GetExpSize(); ++i)
479 exp->GetPointsKeys();
481 exp->GetMetricInfo()->GetJac(pkey);
483 int offset =
m_fields[0]->GetPhys_Offset(i);
485 if (exp->GetMetricInfo()->GetGtype() ==
500 else if (tempEval ==
"Metric")
503 m_graph->GetAllQuadGeoms().size() == 0) ||
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.");
515 for (nv = 0; nv < nVel; ++nv)
521 for (i = 0; i < nVel; ++i)
528 for (i = 0; i <
m_fields[0]->GetExpSize(); ++i)
535 exp->GetMetricInfo()->GetDeriv(pkey);
536 int offset =
m_fields[0]->GetPhys_Offset(i);
545 for (j = 0; j < deriv[0][0].num_elements(); ++j)
547 for (k = 0; k < nVel; ++k)
549 for (l = 0; l < nVel; ++l)
551 jac(l,k) = deriv[k][l][j];
555 jacIdeal = jac * i2rm;
560 char jobvl =
'N', jobvr =
'V';
561 int worklen = 8*nVel, info;
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);
579 for (nv = 0; nv < nVel; ++nv)
581 eval(nv,nv) =
m_beta * (sqrt(eval(nv,nv)) - 1.0);
584 DNekMat beta = evec * eval * evecinv;
586 for (k = 0; k < nVel; ++k)
588 for (l = 0; l < nVel; ++l)
590 m_stress[k][l][offset+j] = beta(k,l);
595 if (deriv[0][0].num_elements() != exp->GetTotPoints())
598 for (k = 0; k < nVel; ++k)
600 for (l = 0; l < nVel; ++l)
603 exp->GetTotPoints(), m_stress[k][l][offset],
604 tmp = m_stress[k][l] + offset, 1);
612 for (i = 0; i < nVel; ++i)
614 for (j = 0; j < nVel; ++j)
616 m_fields[i]->PhysDeriv(j, m_stress[i][j], tmpderiv);
625 else if (tempEval !=
"None")
627 ASSERTL0(
false,
"Unknown temperature form: " + tempEval);
644 for (nv = 0; nv < nVel; ++nv)
648 m_fields[nv]->IProductWRTBase_IterPerExp(forcing[nv], tmp);
652 for (i = 0; i <
m_fields[nv]->GetExpSize(); ++i)
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);
658 for (j = 0; j < nBnd; ++j)
660 forCoeffs[nVel*offset + nv*nBnd + j] =
663 for (j = 0; j < nInt; ++j)
665 forCoeffs[nVel*(offset + nBnd) + nv*nInt + j] =
677 map<int, vector<MultiRegions::ExtraDirDof> > &extraDirDofs =
679 map<int, vector<MultiRegions::ExtraDirDof> >
::iterator it;
681 for (nv = 0; nv < nVel; ++nv)
684 =
m_fields[nv]->GetBndCondExpansions();
687 for (it = extraDirDofs.begin(); it != extraDirDofs.end(); ++it)
689 for (i = 0; i < it->second.size(); ++i)
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>();
704 for (nv = 0; nv < nVel; ++nv)
707 =
m_fields[nv]->GetBndCondExpansions();
711 for (i = 0; i < bndCondExp.num_elements(); ++i)
713 if (
m_fields[nv]->GetBndConditions()[i]->GetBoundaryConditionType()
717 bndCondExp[i]->GetCoeffs();
719 for (j = 0; j < bndCondExp[i]->GetNcoeffs(); ++j)
724 inout[bndMap[bndcnt++]] = sign * bndCoeffs[j];
729 bndcnt += bndCondExp[i]->GetNcoeffs();
756 for (nv = 0; nv < nVel; ++nv)
758 for (i = 0; i <
m_fields[nv]->GetExpSize(); ++i)
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);
764 for (j = 0; j < nBnd; ++j)
767 tmp[nVel*offset + nv*nBnd + j];
769 for (j = 0; j < nInt; ++j)
772 tmp[nVel*(offset + nBnd) + nv*nInt + j];
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;
815 for (i = 0; i < nVel; ++i)
817 for (j = 0; j < nVel; ++j)
820 for (k = 0; k < nB; ++k)
822 for (l = 0; l < nB; ++l)
824 (*A)(k + i*nB, l + j*nB) =
828 for (l = 0; l < nI; ++l)
830 (*B)(k + i*nB, l + j*nI) =
836 for (k = 0; k < nI; ++k)
838 for (l = 0; l < nB; ++l)
840 (*C)(k + i*nI, l + j*nB) =
844 for (l = 0; l < nI; ++l)
846 (*D)(k + i*nI, l + j*nI) =
856 (*A) = (*A) - (*B)*(*C);
886 const int nCoeffs = exp->GetNcoeffs();
887 const int nPhys = exp->GetTotPoints();
891 nCoeffs, nCoeffs, 0.0,
eFULL);
896 for (i = 0; i < nCoeffs; ++i)
901 exp->BwdTrans ( tmp1, tmp2);
902 exp->PhysDeriv (k1, tmp2, tmp3);
903 exp->IProductWRTDerivBase(k2, tmp3, tmp1);
906 nCoeffs, scale, &tmp1[0], 1, &(ret->GetPtr())[0]+i*nCoeffs, 1);
914 std::vector<std::string> &variables)
916 const int nVel =
m_fields[0]->GetCoordim(0);
917 const int nCoeffs =
m_fields[0]->GetNcoeffs();
924 for (
int i = 0; i < nVel; ++i)
928 fieldcoeffs.push_back(tFwd);
930 "ThermStressDiv" + boost::lexical_cast<std::string>(i));
938 for (
int i = 0; i < nVel; ++i)
940 for (
int j = 0; j < nVel; ++j)
944 fieldcoeffs.push_back(tFwd);
947 + boost::lexical_cast<std::string>(i)
948 + boost::lexical_cast<std::string>(j));
#define ASSERTL0(condition, msg)
virtual void v_InitObject()
Set up the linear elasticity system.
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.
std::vector< PointsKey > PointsKeyVector
A base class for describing how to solve specific equations.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
std::vector< std::pair< std::string, std::string > > SummaryList
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
NekDouble m_nu
Poisson ratio.
virtual void v_DoSolve()
Solve elliptic linear elastic system.
NekDouble m_beta
Parameter dictating amount of thermal stress to add.
boost::shared_ptr< ContField2D > ContField2DSharedPtr
std::map< ConstFactorType, NekDouble > ConstFactorMap
DNekScalBlkMatSharedPtr m_BinvD
Matrix of elemental components.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
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.
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Initialisation object for EquationSystem.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
static PreconditionerSharedPtr NullPreconditionerSharedPtr
This class is the abstraction of a global continuous two- dimensional spectral/hp element expansion w...
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Generate summary at runtime.
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
DNekScalBlkMatSharedPtr m_Dinv
Matrix of elemental components.
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
void Neg(int n, T *x, const int incx)
Negate x = -x.
boost::shared_ptr< Expansion > ExpansionSharedPtr
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.
Describe a linear system.
DNekScalBlkMatSharedPtr m_schurCompl
Schur complement boundary-boundary matrix.
EquationSystemFactory & GetEquationSystemFactory()
DNekScalBlkMatSharedPtr m_C
Matrix of elemental components.
NekDouble m_E
Young's modulus.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
DNekMatSharedPtr BuildLaplacianIJMatrix(const int k1, const int k2, const NekDouble scale, LocalRegions::ExpansionSharedPtr exp)
Construct a LaplacianIJ matrix for a given expansion.
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.
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
void BuildMatrixSystem()
Build matrix system for linear elasticity equations.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
boost::shared_ptr< ContField3D > ContField3DSharedPtr
DNekMat MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Geometry is curved or has non-constant factors.
boost::shared_ptr< Geometry > GeometrySharedPtr
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.