28 return (fespace_ptr->
GetMyRank() == 0 && print_level > 0);
33 os <<
"\n<Boundary Info>\n"
34 <<
" Boundary Conditions:\n";
37 os <<
" Boundary " << it.first <<
": ";
58 bool first_print_statement =
true;
62 os <<
" Inhomogeneous Dirichlet defined on ";
65 if (!first_print_statement)
71 first_print_statement =
false;
73 os << it.first <<
"(=" << it.second <<
")";
77 os <<
"<Boundary Info>\n\n";
89 if (boundary.
Find(it.first) == -1)
91 MFEM_ABORT(
" Boundary "
93 <<
" is not defined on the mesh but in Boundary class."
101 std::vector<int> boundary_attributes_keys;
102 for (
int i = 0; i < boundary.
Size(); i++)
106 boundary_attributes_keys.push_back(boundary[i]);
109 if (!boundary_attributes_keys.empty())
112 for (
const auto &it : boundary_attributes_keys)
116 mfem::out <<
") are defined on the mesh but not in the";
117 mfem::out <<
" boundary attributes (Use Neumann).";
125 MFEM_ABORT(
" Periodic boundaries must be defined on the mesh"
126 <<
", not in Boundaries. Exiting...")
140 mfem::out <<
"<Boundary::ComputeBoundaryError>"
163 mfem::out <<
"->Boundary " << i + 1 <<
"\n";
166 mfem::out <<
" Gamma : " << gamma <<
"\n";
167 mfem::out <<
" Average : " << avg <<
"\n";
168 mfem::out <<
" Error : " << error <<
"\n\n";
174 mfem::out <<
"<Boundary::ComputeBoundaryError>" << std::endl;
246 real_t &nrm = loc_vals[0];
247 real_t &avg = loc_vals[1];
248 real_t &error = loc_vals[2];
254 const bool a_is_zero =
alpha == 0.0;
255 const bool b_is_zero =
beta == 0.0;
258 MFEM_ASSERT(fes.
GetVDim() == 1,
"");
265 for (
int i = 0; i < mesh.
GetNBE(); i++)
280 const int int_order = 2 * fe.
GetOrder() + 3;
313 val +=
beta * (shape * loc_dofs);
317 nrm += ip.
weight * face_weight;
320 avg += val * ip.
weight * face_weight;
324 error += (val * val) * ip.
weight * face_weight;
332 real_t glb_nrm = glb_vals[0];
333 real_t glb_avg = glb_vals[1];
334 glb_err = glb_vals[2];
337 if (std::abs(glb_nrm) > 0.0)
345 glb_err = (glb_err >= 0.0) ? sqrt(glb_err) : -sqrt(-glb_err);
356 fespace_ptr_(fespace),
368 mfem::out <<
"<SPDESolver> Initialize Solver .." << std::endl;
377 if (bdr_attributes.Size() > 0)
380 nbc = bdr_attributes.
Max() - bdr_attributes.Min() + 1;
395 dbc_marker_[it.first - 1] = 1;
398 rbc_marker_[it.first - 1] = 1;
414 alpha_ = (nu_ +
dim / 2.0) / 2.0;
415 integer_order_of_exponent_ =
static_cast<int>(std::floor(alpha_));
416 real_t exponent_to_approximate = alpha_ - integer_order_of_exponent_;
419 ComputeRationalCoefficients(exponent_to_approximate);
424 auto diffusion_tensor =
447 if (prolongation_matrix_)
453 mfem::err <<
"<SPDESolver> prolongation matrix is not defined" << std::endl;
455 if (restriction_matrix_)
461 mfem::err <<
"<SPDESolver> restriction matrix is not defined" << std::endl;
468 <<
" [s]" << std::endl;
487 if (integer_order_of_exponent_ > 0)
491 mfem::out <<
"<SPDESolver> Solving PDE (A)^" << integer_order_of_exponent_
492 <<
" u = f" << std::endl;
494 ActivateRepeatedSolve();
495 Solve(
b, helper_gf, 1.0, 1.0, integer_order_of_exponent_);
508 DeactivateRepeatedSolve();
519 for (
int i = 0; i < coeffs_.
Size(); i++)
523 mfem::out <<
"\n<SPDESolver> Solving PDE -Δ u + " << -poles_[i]
524 <<
" u = " << coeffs_[i] <<
" g " << std::endl;
527 Solve(
b, helper_gf, 1.0 - poles_[i], coeffs_[i]);
542 <<
" [s]" << std::endl;
559 MFEM_ABORT(
"Need to call SPDESolver::SetupRandomFieldGenerator(...) first");
565 (*b_wn) *= normalization;
591 const real_t gamma2 = tgamma(nu);
592 return sqrt(pow(2 * M_PI,
dim / 2.0) * det * gamma1 /
593 (gamma2 * pow(nu,
dim / 2.0)));
604 const real_t c1 = cos(e1);
605 const real_t s1 = sin(e1);
606 const real_t c2 = cos(e2);
607 const real_t s2 = sin(e2);
608 const real_t c3 = cos(e3);
609 const real_t s3 = sin(e3);
613 R(0, 0) = c1 * c3 - c2 * s1 * s3;
614 R(0, 1) = -c1 * s3 - c2 * c3 * s1;
616 R(1, 0) = c3 * s1 + c1 * c2 * s3;
617 R(1, 1) = c1 * c2 * c3 - s1 * s3;
625 l(0) = std::pow(l1, 2);
626 l(1) = std::pow(l2, 2);
627 l(2) = std::pow(l3, 2);
628 l *= (1 / (2.0 * nu));
638 const real_t c1 = cos(e1);
639 const real_t s1 = sin(e1);
646 l(0) = std::pow(l1, 2);
647 l(1) = std::pow(l2, 2);
648 l *= (1 / (2.0 * nu));
656 res(0, 0) = std::pow(l1, 2) / (2.0 * nu);
667 if (prolongation_matrix_)
684 restriction_matrix_->
Mult(x, X_);
688 Add(1.0, stiffness_,
alpha, mass_bc_);
689 HypreParMatrix *Ae = Op->EliminateRowsCols(ess_tdof_list_);
690 Op->EliminateBC(*Ae, ess_tdof_list_, X_, B_);
692 for (
int i = 0; i < exponent; i++)
695 SolveLinearSystem(Op);
702 GridFunctionCoefficient gfc(&x);
703 ParLinearForm previous_solution(fespace_ptr_);
704 previous_solution.AddDomainIntegrator(
new DomainLFIntegrator(gfc));
705 previous_solution.Assemble();
707 Op->EliminateBC(*Ae, ess_tdof_list_, X_, B_);
714void SPDESolver::LiftSolution(ParGridFunction &x)
723 mfem::out <<
"\n<SPDESolver> Applying inhomogeneous DBC" << std::endl;
727 ParGridFunction helper_gf(fespace_ptr_);
731 for (
const auto &bc : bc_.dirichlet_coefficients)
735 marker[bc.first - 1] = 1;
736 ConstantCoefficient cc(bc.second);
737 helper_gf.ProjectBdrCoefficient(cc, marker);
741 ParLinearForm
b(fespace_ptr_);
742 ConstantCoefficient zero(0.0);
743 b.AddDomainIntegrator(
new DomainLFIntegrator(zero));
747 Solve(
b, helper_gf, 1.0, 1.0);
756void SPDESolver::UpdateRHS(ParLinearForm &
b)
const
758 if (!repeated_solve_)
763 if (restriction_matrix_)
775void SPDESolver::SolveLinearSystem(
const HypreParMatrix *Op)
777 HypreBoomerAMG prec(*Op);
778 prec.SetPrintLevel(-1);
779 CGSolver cg(fespace_ptr_->
GetComm());
783 cg.SetPreconditioner(prec);
785 cg.SetPrintLevel(std::max(0, print_level_ - 1));
789void SPDESolver::ComputeRationalCoefficients(
real_t exponent)
791 if (abs(exponent) > 1e-12)
795 mfem::out <<
"<SPDESolver> Approximating the fractional exponent "
796 << exponent << std::endl;
802 alpha_ = exponent + integer_order_of_exponent_;
806 integer_order_ =
true;
809 mfem::out <<
"<SPDESolver> Treating integer order PDE." << std::endl;
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
int Find(const T &el) const
Return the first index where 'el' is found; return -1 if not found.
A coefficient that is constant across space and time.
Data type dense matrix using column-major storage.
void Transpose()
(*this) = (*this)^t
real_t InnerProduct(const real_t *x, const real_t *y) const
Compute y^t A x.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
int GetVDim() const
Returns vector dimension.
Abstract class for all finite elements.
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Class for integration point with weight.
Class for an integration rule - an Array of IntegrationPoint.
int GetNPoints() const
Returns the number of the points in the integration rule.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
A matrix coefficient that is constant in space and time.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
int Dimension() const
Dimension of the reference space used within the elements.
int SpaceDimension() const
Dimension of the physical space containing the mesh.
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
int GetNBE() const
Returns number of boundary elements.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Abstract parallel finite element space.
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
const Operator * GetProlongationMatrix() const override
The returned Operator is owned by the FiniteElementSpace.
void GetElementDofs(int i, Array< int > &dofs, DofTransformation &doftrans) const override
The same as GetElementDofs(), but with a user-allocated DofTransformation object. doftrans must be al...
const SparseMatrix * GetRestrictionMatrix() const override
Get the R matrix which restricts a local dof vector to true dof vector.
ParMesh * GetParMesh() const
const FiniteElement * GetFE(int i) const override
Class for parallel grid function.
ParFiniteElementSpace * ParFESpace() const
Class for parallel meshes.
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
void Start()
Start the stopwatch. The elapsed time is not cleared.
void Stop()
Stop the stopwatch.
void SetSize(int s)
Resize the vector to size s.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
void Solve(ParLinearForm &b, ParGridFunction &x)
static DenseMatrix ConstructMatrixCoefficient(real_t l1, real_t l2, real_t l3, real_t e1, real_t e2, real_t e3, real_t nu, int dim)
static real_t ConstructNormalizationCoefficient(real_t nu, real_t l1, real_t l2, real_t l3, int dim)
void GenerateRandomField(ParGridFunction &x)
void SetupRandomFieldGenerator(int seed=0)
SPDESolver(real_t nu, const Boundary &bc, ParFiniteElementSpace *fespace, real_t l1=0.1, real_t l2=0.1, real_t l3=0.1, real_t e1=0.0, real_t e2=0.0, real_t e3=0.0)
void ComputePartialFractionApproximation(real_t &alpha, Array< real_t > &coeffs, Array< real_t > &poles, real_t lmax=1000., real_t tol=1e-10, int npoints=1000, int max_order=100)
real_t IntegrateBC(const ParGridFunction &x, const Array< int > &bdr, real_t alpha, real_t beta, real_t gamma, real_t &glb_err)
bool PrintOutput(const ParFiniteElementSpace *fespace_ptr, int print_level)
void CalcOrtho(const DenseMatrix &J, Vector &n)
void MultADBt(const DenseMatrix &A, const Vector &D, const DenseMatrix &B, DenseMatrix &ADBt)
ADBt = A D B^t, where D is diagonal.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt = A D A^t, where D is diagonal.
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Helper struct to convert a C++ type to an MPI type.
std::map< int, BoundaryType > boundary_attributes
Map to assign homogeneous boundary conditions to defined boundary types.
std::map< int, real_t > dirichlet_coefficients
Coefficient for inhomogeneous Dirichlet boundary conditions.
void PrintInfo(std::ostream &os=mfem::out) const
Print the information specifying the boundary conditions.
void ComputeBoundaryError(const ParGridFunction &solution)
void AddInhomogeneousDirichletBoundaryCondition(int boundary, real_t coefficient)
Add a inhomogeneous Dirichlet boundary condition to the boundary.
void SetRobinCoefficient(real_t coefficient)
Set the robin coefficient for the boundary.
void UpdateIntegrationCoefficients(int i, real_t &alpha, real_t &beta, real_t &gamma)
void VerifyDefinedBoundaries(const Mesh &mesh) const
void AddHomogeneousBoundaryCondition(int boundary, BoundaryType type)
Add a homogeneous boundary condition to the boundary.