- Note
- Prerequisites of this tutorial include MSH-1: Create a 2D mesh from Gmsh
- Note
- Intended learning outcome:
- general structure of a program developed using MoFEM
- idea of Simple Interface in MoFEM and how to use it
- solving eigen problem
- Use of form integrators
- How to push the developed UDOs to the Pipeline
Introduction
This tutorial presents the solution of an eigen-elastic problem using the operator OpGradSymTensorGrad, together with the SLEPc EPS solver.
Notation description
Contraction: given a fourth-order tensor \( A_{n_1 \times n_2 \times n_3 \times n_4} \) and a matrix \( B_{n_3 \times n_4} \), their contraction is written as \( C = A : B \), with
\[
C_{ij} = \sum_{k l} A_{i j k l} B_{k l}.
\label{contraction}
\]
Einstein notation:
repeated indices imply summation. For instance,
\[
\label{einstein}
y=\sum^{n}_{i=1}a_{i}x_{i}\quad \Leftrightarrow \quad y=a_ix_i
\]
Using this rule, the contraction above becomes
\[
C_{ij}=A_{ijkl}B_{kl}
\]
Problem statement
Strong form
For a time-dependent elastic problem, the displacement field \(\mathbf{u}\) and the stress \(\pmb{\sigma}\) satisfy balance of linear momentum:
\[
\label{eq:time_elastic}
\rho \frac{\partial^{2} \mathbf{u}}{\partial t^{2}} = \nabla \cdot \pmb{\sigma}\quad \text{in } \Omega
\]
where \(\rho\) is the material density.
The stress is related to strain via tensor contraction:
\[
\label{eq:stress-strain}
{\pmb{\sigma}(\mathbf{u})}=\mathbf{C}:\pmb{\varepsilon}(\mathbf{u})
\]
with the strain
\[
\begin{equation}
\varepsilon(\mathbf{u}) = \frac{1}{2} \big(\nabla \mathbf{u} + \nabla \mathbf{u}^T \big).
\end{equation}
\]
The stiffness tensor \(\mathbf{C}\) depends on Young's modulus \(E\) and Poisson ratio \(\nu\).
As this problem focuses on determining the natural frequencies of a tuning fork in free vibration, which is treated as an eigenvalue problem later, the displacement \(\mathbf{u}(\mathbf{x},t)\) is decomposed into position-dependent part and time-dependent phase:
\[
\label{eq:time_harmonic}
\mathbf{u}(\mathbf{x},t) = \mathbf{\hat{u}}(\mathbf{x}) e^{i \omega t}
\]
Combining Eqs. \eqref{eq:time_elastic} and \eqref{eq:time_harmonic} leads to:
\[
\label{eq:static}
\nabla \cdot \pmb{\sigma}(\mathbf{\hat{u}}) + \rho \omega^2 \mathbf{\hat{u}} = 0
\]
The bulk modulus \(K\) and the shear modulus \(G\) are related to \(E\) and \(\nu\) by:
\[
\label{eq:bulk_mod}
K = \frac{E}{3(1-2\nu)}, \quad
G = \frac{E}{2(1+\nu)}
\]
Then the stiffness tensor is
\[
\label{eq:stiffness_coe}
\mathbf{C}_{ijkl} = 2G \delta_{ik}\delta_{jl} + A \Big(K - \frac{2}{3} G \Big) \delta_{ij} \delta_{kl}
\]
where \(A = \frac{2G}{K+\frac{4}{3}G}\) for 2D plane strain and \(A = 1\) for 3D.
Therefore we get the final strong form of the problem:
\[
\begin{align}
\label{eq:strong_form_domain}
\nabla \!\cdot\! \left[\mathbf{C} : \varepsilon(\hat{\mathbf{u}})\right] + \rho \omega^2 \hat{\mathbf{u}} &= 0, && \text{in } \Omega, \\[2mm]
\label{eq:strong_form_boundary}
\pmb{\sigma}(\hat{\mathbf{u}}) \cdot \mathbf{n}&= 0, && \text{on } \Gamma.
\end{align}
\]
Weak form
Denote \(\delta\mathbf{\hat{u}}\) as the test function for the displacement. As \(\mathbf{C}:\varepsilon(\hat{\mathbf{u}})\) is symmetric, we have
\[
\nabla \delta \mathbf{\hat{u}}:\mathbf{C}:\varepsilon(\hat{\mathbf{u}})=\varepsilon(\delta \mathbf{\hat{u}}) : \mathbf{C} : \varepsilon(\mathbf{\hat{u}})
\]
Starting from Eq. \eqref{eq:strong_form_domain} and \eqref{eq:strong_form_boundary}, the weak form is:
Find \(\mathbf{\hat{u}} \in \mathrm{H}^1(\Omega)\) such that
\[
\label{eq:weak_form}
\int_{\Omega} \varepsilon(\delta \mathbf{\hat{u}}) : \mathbf{C} : \varepsilon(\mathbf{\hat{u}}) \, d\Omega -\int_{\Gamma}
(\delta \mathbf{\hat{u}})\cdot \left[\mathbf{C} : \varepsilon(\mathbf{\hat{u}})\right] \cdot \mathbf{n} \, d\Gamma
= \omega^2 \int_{\Omega} \rho (\delta \mathbf{\hat{u}})\cdot \mathbf{\hat{u}} \, d\Omega,
\quad \forall \delta \mathbf{\hat{u}} \in \mathrm{H}^1(\Omega)
\]
According to the natural boundary condition in Eq. \eqref{eq:strong_form_boundary}, the second term on the left-hand side vanishes.
Approximating test and trial functions using shape functions \(\pmb{\phi}\)
\[
\label{eq:shape_functions}
\begin{align}
\delta\mathbf{\hat{u}}\approx\delta\mathbf{\hat{u}}^h&=\sum_{i=1}^{n}\pmb{\phi}_i\delta\mathbf{\hat{u}}_i, \\[4pt]
\mathbf{\hat{u}}\approx\mathbf{\hat{u}}^h&=\sum_{j=1}^{n}\pmb{\phi}_j\mathbf{\hat{u}}_j
\end{align}
\]
where \(n\) is the total number of DOFs, we obtain the discrete form:
\[
\label{eq:discre_weak_form}
\sum_{j=1}^n \int_{\Omega} \varepsilon(\pmb{\phi}_i) : \mathbf{C} : \varepsilon(\pmb{\phi}_j) \, d\Omega \, \mathbf{\hat{u}}_j
= \sum_{j=1}^n \omega^2 \int_{\Omega} \rho \pmb{\phi}_i \pmb{\phi}_j \, d\Omega \, \mathbf{\hat{u}}_j, \quad \forall i = 1,2,...,n
\]
The stiffness and mass matrices are defined as:
\[
\label{eq:K_matrix}
\mathbf{K}_{ij} = \int_{\Omega} \varepsilon(\pmb{\phi}_i) : \mathbf{C} : \varepsilon(\pmb{\phi}_j) \, d\Omega
\]
\[
\label{eq:M_matrix}
\mathbf{M}_{ij} = \int_{\Omega} \rho \pmb{\phi}_i^T \pmb{\phi}_j \, d\Omega
\]
To distinguish the material parameter from the stiffness matrix, \(\textbf{K}\) is adopted for the stiffness matrix.
This leads to the generalized eigenvalue problem:
\[
\begin{equation}
\mathbf{K} \mathbf{\hat{u}} = \omega^2 \mathbf{M} \mathbf{\hat{u}}
\label{eq:generalized_eigen}
\end{equation}
\]
where \(\lambda = \omega^2\) are the eigenvalues.
Rigid body motion modes
\[
\label{eq:rigid_motion}
\begin{aligned}
\mathbf{r_1} &=
\begin{bmatrix} 1\\0\\0\\1\\0\\0\\\vdots \end{bmatrix}_{3n\times1}, \quad
\mathbf{r_2} =
\begin{bmatrix} 0\\1\\0\\0\\1\\0\\\vdots \end{bmatrix}_{3n\times1}, \quad
\mathbf{r_3} =
\begin{bmatrix} 0\\0\\1\\0\\0\\1\\\vdots \end{bmatrix}_{3n\times1}, \\[2mm]
\mathbf{r_4} &=
\begin{bmatrix} 0\\-z_1\\y_1\\0\\-z_2\\y_2\\\vdots \end{bmatrix}_{3n\times1}, \quad
\mathbf{r_5} =
\begin{bmatrix} z_1\\0\\-x_1\\z_2\\0\\-x_2\\\vdots \end{bmatrix}_{3n\times1}, \quad
\mathbf{r_6} =
\begin{bmatrix} -y_1\\x_1\\0\\-y_2\\x_2\\0\\\vdots \end{bmatrix}_{3n\times1}.
\end{aligned}
\]
In free vibration problems, rigid body modes naturally appear with zero frequency. These modes satisfy \(\mathbf{K}\mathbf{r}=0\), indicating that \(\mathbf{r}\) represent the eigenvectors of \(\mathbf{K}\) with eigenvalues of 0, which can be derived intuitively without solving the eigen problem shown in Eq. \eqref{eq:generalized_eigen}. Above expressions of \(\mathbf{r_1},\dots,\mathbf{r_6}\) only exhibit the values for vertices in \(\mathbf{\hat{u}}\). In a higher-order case, the values for higher-order entities are 0, which are placed after the vertex values. Here \(n\) is the number of mesh vertices. Each vertex has 3 DOFs, so vectors are with the size of \(3n\times1\). The vectors \(\mathbf{r_1}\), \(\mathbf{r_2}\) and \(\mathbf{r_3}\) stand for translation modes, and \(\mathbf{r_4}\), \(\mathbf{r_5}\) and \(\mathbf{r_6}\) stand for rotation modes.
Implementation
Definition of operators and physical parameters
Since this is a generalized eigenvalue problem, we use the EPS solver from SLEPc rather than KSP. ElementsAndOps defines a class DomainEle that later defines the operators required regardless of the space dimension:
OpDomainGradGrad assembles the stiffness matrix \(\mathbf{K}\) (see Eq. \eqref{eq:K_matrix}), while OpDomainMass assembles the mass matrix \(\mathbf{M}\) (see Eq. \eqref{eq:M_matrix}).
};
using DomainEle = VolumeElementForcesAndSourcesCore;
};
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
#define EXECUTABLE_DIMENSION
double young_modulus
Young modulus.
double poisson_ratio
Poisson ratio.
Example structure and problem setup
The major functions required for the FEM implementation are readMesh(), setupProblem(), createCommonData(), boundaryCondition(), assembleSystem(), solveSystem(), outputResults() and checkResults(). The matrices \(\mathbf{M}\) and \(\mathbf{K}\) correspond to those defined in Eq. \eqref{eq:K_matrix} and Eq. \eqref{eq:M_matrix}, respectively. The pointer matDPtr is used to construct the fourth-order tensor \(\mathbf{C}\) in Eq. \eqref{eq:stiffness_coe}. rigidBodyMotion stores the six eigenvectors associated with the rigid body motion modes described in Eq. \eqref{eq:rigid_motion}. Each self-defined function in MoFEM starts with MoFEMFunctionBegin and ends with MoFEMFunctionReturn(0), with MoFEMErrorCode serving as the standard return type in MoFEM to report the success or failure of the function.
private:
boost::shared_ptr<MatrixDouble>
matDPtr;
};
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Set up problem]
boost::shared_ptr< MatrixDouble > matDPtr
std::array< SmartPetscObj< Vec >, 6 > rigidBodyMotion
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
Reference to MoFEM interface.
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode outputResults()
[Solve]
Deprecated interface functions.
The Simple interface is called to load the mesh data, and register the solution field (displacement \(\mathbf{\hat{u}^h}\)) with the polynomial base AINSWORTH_BERNSTEIN_BEZIER_BASE approximated in the \(\mathrm{H}^1\) space.
}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
}
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Create common data
The default values of \(\rho\), \(E\) and \(\nu\) can be overriden from either the param_file.petsc file or the command line (if provided) by users, before the dependent variables \(K\) and \(G\) are calculated. Then the stiffness tensor \(\mathbf{C}\) is created using FTensor in the code. The relevant expressions are given in Eq. \eqref{eq:bulk_mod} - \eqref{eq:stiffness_coe}.
auto set_matrial_stiffens = [&]() {
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
matDPtr);
: 1;
};
matDPtr = boost::make_shared<MatrixDouble>();
CHKERR set_matrial_stiffens();
}
Kronecker Delta class symmetric.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Apply boundary conditions
The six rigid body motion vectors are constructed as explained in Eq. \eqref{eq:rigid_motion}, with each assigned a value to every vertex DOF. getDofCoeffIdx() returns the coordinate index \(x,y\) or \(z\), while getPetscGlobalDofIdx() retrieves the corresponding global DOF index. lo_uid is used to access all the vertices in the mesh. The subsequent loop synchronizes the parallel processes by VecAssemblyBegin and VecAssemblyEnd. The functions VecGhostUpdateBegin and VecGhostUpdateEnd handle the exchange of information on the interfaces where data is not shared between adjacent entities.
for (
int n = 1;
n != 6; ++
n)
auto lo_uid =
DofEntity::getUniqueIdCalculate(0, get_id_for_min_type<MBVERTEX>());
auto hi = dofs->upper_bound(lo_uid);
std::array<double, 3> coords;
for (auto lo = dofs->lower_bound(lo_uid); lo != hi; ++lo) {
auto ent = (*lo)->getEnt();
if ((*lo)->getDofCoeffIdx() == 0) {
INSERT_VALUES);
-coords[1], INSERT_VALUES);
-coords[2], INSERT_VALUES);
} else if ((*lo)->getDofCoeffIdx() == 1) {
INSERT_VALUES);
coords[0], INSERT_VALUES);
-coords[2], INSERT_VALUES);
} else if ((*lo)->getDofCoeffIdx() == 2) {
1, INSERT_VALUES);
coords[0], INSERT_VALUES);
coords[1], INSERT_VALUES);
}
}
}
}
SCATTER_FORWARD);
SCATTER_FORWARD);
}
}
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
const double n
refractive index of diffusive medium
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
virtual moab::Interface & get_moab()=0
virtual int get_comm_rank() const =0
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Assemble system
AddHOOps undertakes the same set of tasks as the following combination of individual operations:
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetHOInvJacToScalarBases<SPACE_DIM>(
H1, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetHOWeights(det_ptr));
The calculation above involves the Jacobian matrix \(J\) which maps the gradients expressed in the parent coordinate system to those in the global coordinate system, which are included in AddHOOps now. In this case, the two fields \(\delta\mathbf{\hat{u}}\) and \(\mathbf{\hat{u}}\) passed into OpDomainGradGrad below correspond to the shape functions \(\pmb{\phi}_i\) and \(\pmb{\phi}_j\) in Eq. \eqref{eq:K_matrix}. As \(\varepsilon\) serves as a function mapping the fields to their gradients, OpDomainGradGrad is doing exactly what is shown in Eq. \eqref{eq:K_matrix}. Similarly, OpDomainMass performs the operation described in Eq. \eqref{eq:M_matrix}. The difference in the integration order between the two operators arises from the presence of the two gradient operations in OpDomainGradGrad. Its quadrature order is therefore reduced by two relative to that of OpDomainMass. As we are dealing with a generalized eigenvalue problem, \(M\) and \(K\) need to be explicitly assigned into PipelineManager by pipeline_mng->getDomainLhsFE()->B = M and pipeline_mng->getDomainLhsFE()->B = K. MatZeroEntries is called to clear the matrix before assembly, and MatAssemblyBegin and MatAssemblyEnd are used to finalize the assembly process.
auto calculate_stiffness_matrix = [&]() {
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
};
pipeline_mng->getDomainLhsFE()->B =
K;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR MatAssemblyBegin(
K, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
K, MAT_FINAL_ASSEMBLY);
};
auto calculate_mass_matrix = [&]() {
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
};
pipeline_mng->getDomainLhsFE()->B =
M;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR MatAssemblyBegin(
M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
M, MAT_FINAL_ASSEMBLY);
};
CHKERR calculate_stiffness_matrix();
CHKERR calculate_mass_matrix();
}
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
static constexpr int approx_order
Solve system
An EPS solver is applied with proper configurations applied. EPSSetProblemType(ePS, EPS_GHEP) specifies the problem is a generalized Hermitian eigenvalue problem, and EPS_SMALLEST_MAGNITUDE in EPSSetWhichEigenpairs tells only the eigenpairs with the smallest eigenvalues are required for computing. The six rigid body motion modes are eliminated by declaring them as a nullspace by EPSSetDeflationSpace. The print_info uses MOFEM_LOG_C for desired terminal output. Using plain std::cout is discouraged because it will not handle output properly when running in parallel.
auto create_eps = [](MPI_Comm comm) {
return SmartPetscObj<EPS>(
eps);
};
auto deflate_vectors = [&]() {
std::array<Vec, 6> deflate_vectors;
for (
int n = 0;
n != 6; ++
n) {
}
CHKERR EPSSetDeflationSpace(
ePS, 6, &deflate_vectors[0]);
};
auto print_info = [&]() {
ST st;
PetscInt nev, maxit, its;
" Number of iterations of the method: %d", its);
CHKERR EPSGetDimensions(
ePS, &nev, NULL, NULL);
MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
" Number of requested eigenvalues: %d",
nev);
" Stopping condition: tol=%.4g, maxit=%d", (
double)
tol, maxit);
PetscScalar eigr, eigi;
for (int nn = 0; nn < nev; nn++) {
CHKERR EPSGetEigenpair(
ePS, nn, &eigr, &eigi, PETSC_NULLPTR,
PETSC_NULLPTR);
" ncov = %d eigr = %.4g eigi = %.4g (inv eigr = %.4g)", nn,
eigr, eigi, 1. / eigr);
}
};
auto setup_eps = [&]() {
CHKERR EPSSetWhichEigenpairs(
ePS, EPS_SMALLEST_MAGNITUDE);
};
}
#define MOFEM_LOG_C(channel, severity, format,...)
virtual MPI_Comm & get_comm() const =0
Postprocess
This part handles postprocessing and writing the results to output files. Three quantities are written, displacement \(u\), strain \(\varepsilon\) and stress \(\sigma\). These are accessed through pointers such as u_ptr, strain_ptr, and stress_ptr. Statements like post_proc_fe->getOpPtrVector().push_back(new Opxxx(....)) add postprocessing operators that take the FEM solution and compute derived quantities, such as strain and stress. The operators Opxxxx() can either be built-in or user-defined. The call post_proc_fe->getOpPtrVector().push_back(new OpPPMap()) maps the output variables onto the mesh. The function EPSGetEigenpair retrives the eigenvalues and eigenvectors and stores them in eigr (real part of the eigenvalues), eigi (imaginary part of the eigenvalues) and D (real part of the eigenvectors). The imaginary part of eigenvectors are omitted as displacement in the real world can't be complex numbers. DMoFEMMeshToLocalVector transfers data either from a vector to the mesh or from the mesh to a vector, depending on whether the flag is SCATTER_REVERSE or SCATTER_FORWARD.
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
post_proc_fe->getOpPtrVector(), {H1});
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
auto strain_ptr = boost::make_shared<MatrixDouble>();
auto stress_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U", grad_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpSymmetrizeTensor<SPACE_DIM>(grad_ptr, strain_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{},
OpPPMap::DataMapMat{{"U", u_ptr}},
)
);
pipeline_mng->getDomainPostProcFE() = post_proc_fe;
PetscInt nev;
CHKERR EPSGetDimensions(ePS, &nev, NULL, NULL);
PetscScalar eigr, eigi, nrm2r;
for (int nn = 0; nn < nev; nn++) {
CHKERR EPSGetEigenpair(ePS, nn, &eigr, &eigi,
D, PETSC_NULLPTR);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
" ncov = %d omega2 = %.8g omega = %.8g frequency = %.8g", nn,
eigr, std::sqrt(std::abs(eigr)),
std::sqrt(std::abs(eigr)) / (2 * M_PI));
CHKERR pipeline_mng->loopFiniteElementsPostProc();
post_proc_fe->writeFile("out_eig_" + boost::lexical_cast<std::string>(nn) +
".h5m");
}
}
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Validate results
Check if the first eigenvalue matches the analytically predicted value.
PetscBool test_flg = PETSC_FALSE;
PETSC_NULLPTR);
if (test_flg) {
PetscScalar eigr, eigi;
CHKERR EPSGetEigenpair(
ePS, 0, &eigr, &eigi, PETSC_NULLPTR, PETSC_NULLPTR);
constexpr double regression_value = 12579658;
if (fabs(eigr - regression_value) > 1)
"Regression test faileed; wrong eigen value.");
}
}
@ MOFEM_ATOM_TEST_INVALID
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Main function
This is the main function of the tutorial. No modifications are needed.
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
SlepcInitialize(&argc, &argv, param_file,
help);
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
CHKERR DMRegister_MoFEM(dm_name);
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
SlepcFinalize();
}
#define CATCH_ERRORS
Catch errors.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Results
The mesh is shown in Figure 1. The vibration modes corresponding to the six smallest natural frequencies are presented in Figure 2(a)-(f). A comparison between the Calculix and our implementation is given in Table 1. As shown, the deviations in our results are less than 1%, demonstrating good agreement between the two approaches.
Figure 1: Geometry of the problem
(a) Mode 1 |
(b) Mode 2 |
(c) Mode 3 |
(d) Mode 4 |
(e) Mode 5 |
(f) Mode 6 |
Figure 2: The 1st to 6th vibration modes for the tuning fork.
Table 1: Natural frequencies (Hz) of the tuning fork modes
| Mode | Calculix | MoFEM |
| 1 | 440.05 | 442.98 |
| 2 | 673.51 | 676.92 |
| 3 | 1689.30 | 1692.00 |
| 4 | 1825.55 | 1826.71 |
| 5 | 2777.73 | 2790.89 |
Source code
#undef EPS
#include <slepceps.h>
};
};
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
private:
boost::shared_ptr<MatrixDouble>
matDPtr;
};
auto set_matrial_stiffens = [&]() {
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
matDPtr);
: 1;
};
matDPtr = boost::make_shared<MatrixDouble>();
CHKERR set_matrial_stiffens();
}
}
}
}
for (
int n = 1;
n != 6; ++
n)
auto lo_uid =
DofEntity::getUniqueIdCalculate(0, get_id_for_min_type<MBVERTEX>());
auto hi = dofs->upper_bound(lo_uid);
std::array<double, 3> coords;
for (auto lo = dofs->lower_bound(lo_uid); lo != hi; ++lo) {
auto ent = (*lo)->getEnt();
if ((*lo)->getDofCoeffIdx() == 0) {
INSERT_VALUES);
-coords[1], INSERT_VALUES);
-coords[2], INSERT_VALUES);
} else if ((*lo)->getDofCoeffIdx() == 1) {
INSERT_VALUES);
coords[0], INSERT_VALUES);
-coords[2], INSERT_VALUES);
} else if ((*lo)->getDofCoeffIdx() == 2) {
1, INSERT_VALUES);
coords[0], INSERT_VALUES);
coords[1], INSERT_VALUES);
}
}
}
}
SCATTER_FORWARD);
SCATTER_FORWARD);
}
}
auto calculate_stiffness_matrix = [&]() {
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
};
pipeline_mng->getDomainLhsFE()->B =
K;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR MatAssemblyBegin(
K, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
K, MAT_FINAL_ASSEMBLY);
};
auto calculate_mass_matrix = [&]() {
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
};
pipeline_mng->getDomainLhsFE()->B =
M;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR MatAssemblyBegin(
M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
M, MAT_FINAL_ASSEMBLY);
};
CHKERR calculate_stiffness_matrix();
CHKERR calculate_mass_matrix();
}
auto create_eps = [](MPI_Comm comm) {
};
auto deflate_vectors = [&]() {
std::array<Vec, 6> deflate_vectors;
for (
int n = 0;
n != 6; ++
n) {
}
CHKERR EPSSetDeflationSpace(
ePS, 6, &deflate_vectors[0]);
};
auto print_info = [&]() {
ST st;
PetscInt nev, maxit, its;
" Number of iterations of the method: %d", its);
CHKERR EPSGetDimensions(
ePS, &nev, NULL, NULL);
MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
" Number of requested eigenvalues: %d",
nev);
" Stopping condition: tol=%.4g, maxit=%d", (
double)
tol, maxit);
PetscScalar eigr, eigi;
for (int nn = 0; nn < nev; nn++) {
CHKERR EPSGetEigenpair(
ePS, nn, &eigr, &eigi, PETSC_NULLPTR,
PETSC_NULLPTR);
" ncov = %d eigr = %.4g eigi = %.4g (inv eigr = %.4g)", nn,
eigr, eigi, 1. / eigr);
}
};
auto setup_eps = [&]() {
CHKERR EPSSetWhichEigenpairs(
ePS, EPS_SMALLEST_MAGNITUDE);
};
}
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
post_proc_fe->getOpPtrVector(), {H1});
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
auto strain_ptr = boost::make_shared<MatrixDouble>();
auto stress_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{},
OpPPMap::DataMapMat{{"U", u_ptr}},
)
);
pipeline_mng->getDomainPostProcFE() = post_proc_fe;
PetscInt nev;
CHKERR EPSGetDimensions(ePS, &nev, NULL, NULL);
PetscScalar eigr, eigi, nrm2r;
for (int nn = 0; nn < nev; nn++) {
CHKERR EPSGetEigenpair(ePS, nn, &eigr, &eigi,
D, PETSC_NULLPTR);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
" ncov = %d omega2 = %.8g omega = %.8g frequency = %.8g", nn,
eigr, std::sqrt(std::abs(eigr)),
std::sqrt(std::abs(eigr)) / (2 * M_PI));
CHKERR pipeline_mng->loopFiniteElementsPostProc();
post_proc_fe->writeFile("out_eig_" + boost::lexical_cast<std::string>(nn) +
".h5m");
}
}
PetscBool test_flg = PETSC_FALSE;
PETSC_NULLPTR);
if (test_flg) {
PetscScalar eigr, eigi;
CHKERR EPSGetEigenpair(
ePS, 0, &eigr, &eigi, PETSC_NULLPTR, PETSC_NULLPTR);
constexpr double regression_value = 12579658;
if (fabs(eigr - regression_value) > 1)
"Regression test faileed; wrong eigen value.");
}
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
SlepcInitialize(&argc, &argv, param_file,
help);
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
SlepcFinalize();
}
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
implementation of Data Operators for Forces and Sources
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for MatrixDouble vector field values calculation.
Post post-proc data at points from hash maps.
Operator for symmetrizing tensor fields.
PipelineManager interface.
Simple interface for fast problem set-up.
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode getOptions()
get options
intrusive_ptr for managing petsc objects
Volume finite element base.
