909 lines
22 KiB
C++
909 lines
22 KiB
C++
// MFEM Example 37 - Serial/Parallel Shared Code
|
||
|
||
#include "mfem.hpp"
|
||
#include <fstream>
|
||
#include <iostream>
|
||
#include <functional>
|
||
|
||
namespace mfem
|
||
{
|
||
|
||
/// @brief Inverse sigmoid function
|
||
real_t inv_sigmoid(real_t x)
|
||
{
|
||
real_t tol = 1e-12;
|
||
x = std::min(std::max(tol,x), real_t(1.0)-tol);
|
||
return std::log(x/(1.0-x));
|
||
}
|
||
|
||
/// @brief Sigmoid function
|
||
real_t sigmoid(real_t x)
|
||
{
|
||
if (x >= 0)
|
||
{
|
||
return 1.0/(1.0+std::exp(-x));
|
||
}
|
||
else
|
||
{
|
||
return std::exp(x)/(1.0+std::exp(x));
|
||
}
|
||
}
|
||
|
||
/// @brief Derivative of sigmoid function
|
||
real_t der_sigmoid(real_t x)
|
||
{
|
||
real_t tmp = sigmoid(-x);
|
||
return tmp - std::pow(tmp,2);
|
||
}
|
||
|
||
/// @brief Returns f(u(x)) where u is a scalar GridFunction and f:R → R
|
||
class MappedGridFunctionCoefficient : public GridFunctionCoefficient
|
||
{
|
||
protected:
|
||
std::function<real_t(const real_t)> fun; // f:R → R
|
||
public:
|
||
MappedGridFunctionCoefficient()
|
||
:GridFunctionCoefficient(),
|
||
fun([](real_t x) {return x;}) {}
|
||
MappedGridFunctionCoefficient(const GridFunction *gf,
|
||
std::function<real_t(const real_t)> fun_,
|
||
int comp=1)
|
||
:GridFunctionCoefficient(gf, comp),
|
||
fun(fun_) {}
|
||
|
||
|
||
real_t Eval(ElementTransformation &T,
|
||
const IntegrationPoint &ip) override
|
||
{
|
||
return fun(GridFunctionCoefficient::Eval(T, ip));
|
||
}
|
||
void SetFunction(std::function<real_t(const real_t)> fun_) { fun = fun_; }
|
||
};
|
||
|
||
|
||
/// @brief Returns f(u(x)) - f(v(x)) where u, v are scalar GridFunctions and f:R → R
|
||
class DiffMappedGridFunctionCoefficient : public GridFunctionCoefficient
|
||
{
|
||
protected:
|
||
const GridFunction *OtherGridF;
|
||
GridFunctionCoefficient OtherGridF_cf;
|
||
std::function<real_t(const real_t)> fun; // f:R → R
|
||
public:
|
||
DiffMappedGridFunctionCoefficient()
|
||
:GridFunctionCoefficient(),
|
||
OtherGridF(nullptr),
|
||
OtherGridF_cf(),
|
||
fun([](real_t x) {return x;}) {}
|
||
DiffMappedGridFunctionCoefficient(const GridFunction *gf,
|
||
const GridFunction *other_gf,
|
||
std::function<real_t(const real_t)> fun_,
|
||
int comp=1)
|
||
:GridFunctionCoefficient(gf, comp),
|
||
OtherGridF(other_gf),
|
||
OtherGridF_cf(OtherGridF),
|
||
fun(fun_) {}
|
||
|
||
real_t Eval(ElementTransformation &T,
|
||
const IntegrationPoint &ip) override
|
||
{
|
||
const real_t value1 = fun(GridFunctionCoefficient::Eval(T, ip));
|
||
const real_t value2 = fun(OtherGridF_cf.Eval(T, ip));
|
||
return value1 - value2;
|
||
}
|
||
void SetFunction(std::function<real_t(const real_t)> fun_) { fun = fun_; }
|
||
};
|
||
|
||
/// @brief Solid isotropic material penalization (SIMP) coefficient
|
||
class SIMPInterpolationCoefficient : public Coefficient
|
||
{
|
||
protected:
|
||
GridFunction *rho_filter;
|
||
real_t min_val;
|
||
real_t max_val;
|
||
real_t exponent;
|
||
|
||
public:
|
||
SIMPInterpolationCoefficient(GridFunction *rho_filter_, real_t min_val_= 1e-6,
|
||
real_t max_val_ = 1.0, real_t exponent_ = 3)
|
||
: rho_filter(rho_filter_), min_val(min_val_), max_val(max_val_),
|
||
exponent(exponent_) { }
|
||
|
||
real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override
|
||
{
|
||
real_t val = rho_filter->GetValue(T, ip);
|
||
real_t coeff = min_val + pow(val,exponent)*(max_val-min_val);
|
||
return coeff;
|
||
}
|
||
};
|
||
|
||
|
||
/// @brief Strain energy density coefficient
|
||
class StrainEnergyDensityCoefficient : public Coefficient
|
||
{
|
||
protected:
|
||
Coefficient * lambda=nullptr;
|
||
Coefficient * mu=nullptr;
|
||
GridFunction *u = nullptr; // displacement
|
||
GridFunction *rho_filter = nullptr; // filter density
|
||
DenseMatrix grad; // auxiliary matrix, used in Eval
|
||
real_t exponent;
|
||
real_t rho_min;
|
||
|
||
public:
|
||
StrainEnergyDensityCoefficient(Coefficient *lambda_, Coefficient *mu_,
|
||
GridFunction * u_, GridFunction * rho_filter_, real_t rho_min_=1e-6,
|
||
real_t exponent_ = 3.0)
|
||
: lambda(lambda_), mu(mu_), u(u_), rho_filter(rho_filter_),
|
||
exponent(exponent_), rho_min(rho_min_)
|
||
{
|
||
MFEM_ASSERT(rho_min_ >= 0.0, "rho_min must be >= 0");
|
||
MFEM_ASSERT(rho_min_ < 1.0, "rho_min must be < 1");
|
||
MFEM_ASSERT(u, "displacement field is not set");
|
||
MFEM_ASSERT(rho_filter, "density field is not set");
|
||
}
|
||
|
||
real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override
|
||
{
|
||
real_t L = lambda->Eval(T, ip);
|
||
real_t M = mu->Eval(T, ip);
|
||
u->GetVectorGradient(T, grad);
|
||
real_t div_u = grad.Trace();
|
||
real_t density = L*div_u*div_u;
|
||
int dim = T.GetSpaceDim();
|
||
for (int i=0; i<dim; i++)
|
||
{
|
||
for (int j=0; j<dim; j++)
|
||
{
|
||
density += M*grad(i,j)*(grad(i,j)+grad(j,i));
|
||
}
|
||
}
|
||
real_t val = rho_filter->GetValue(T,ip);
|
||
|
||
return -exponent * pow(val, exponent-1.0) * (1-rho_min) * density;
|
||
}
|
||
};
|
||
|
||
/// @brief Volumetric force for linear elasticity
|
||
class VolumeForceCoefficient : public VectorCoefficient
|
||
{
|
||
private:
|
||
real_t r;
|
||
Vector center;
|
||
Vector force;
|
||
public:
|
||
VolumeForceCoefficient(real_t r_,Vector & center_, Vector & force_) :
|
||
VectorCoefficient(center_.Size()), r(r_), center(center_), force(force_) { }
|
||
|
||
using VectorCoefficient::Eval;
|
||
|
||
void Eval(Vector &V, ElementTransformation &T,
|
||
const IntegrationPoint &ip) override
|
||
{
|
||
Vector xx; xx.SetSize(T.GetDimension());
|
||
T.Transform(ip,xx);
|
||
for (int i=0; i<xx.Size(); i++)
|
||
{
|
||
xx[i]=xx[i]-center[i];
|
||
}
|
||
|
||
real_t cr=xx.Norml2();
|
||
V.SetSize(T.GetDimension());
|
||
if (cr <= r)
|
||
{
|
||
V = force;
|
||
}
|
||
else
|
||
{
|
||
V = 0.0;
|
||
}
|
||
}
|
||
|
||
void Set(real_t r_,Vector & center_, Vector & force_)
|
||
{
|
||
r=r_;
|
||
center = center_;
|
||
force = force_;
|
||
}
|
||
};
|
||
|
||
/**
|
||
* @brief Class for solving Poisson's equation:
|
||
*
|
||
* - ∇ ⋅(κ ∇ u) = f in Ω
|
||
*
|
||
*/
|
||
class DiffusionSolver
|
||
{
|
||
private:
|
||
Mesh * mesh = nullptr;
|
||
int order = 1;
|
||
// diffusion coefficient
|
||
Coefficient * diffcf = nullptr;
|
||
// mass coefficient
|
||
Coefficient * masscf = nullptr;
|
||
Coefficient * rhscf = nullptr;
|
||
Coefficient * essbdr_cf = nullptr;
|
||
Coefficient * neumann_cf = nullptr;
|
||
VectorCoefficient * gradient_cf = nullptr;
|
||
|
||
// FEM solver
|
||
int dim;
|
||
FiniteElementCollection * fec = nullptr;
|
||
FiniteElementSpace * fes = nullptr;
|
||
Array<int> ess_bdr;
|
||
Array<int> ess_tdof_list;
|
||
Array<int> neumann_bdr;
|
||
GridFunction * u = nullptr;
|
||
LinearForm * b = nullptr;
|
||
BilinearForm * a = nullptr;
|
||
OperatorPtr A;
|
||
bool parallel;
|
||
#ifdef MFEM_USE_MPI
|
||
ParMesh * pmesh = nullptr;
|
||
ParFiniteElementSpace * pfes = nullptr;
|
||
#endif
|
||
|
||
public:
|
||
DiffusionSolver() { }
|
||
DiffusionSolver(Mesh * mesh_, int order_, Coefficient * diffcf_,
|
||
Coefficient * cf_);
|
||
|
||
void SetMesh(Mesh * mesh_)
|
||
{
|
||
mesh = mesh_;
|
||
parallel = false;
|
||
#ifdef MFEM_USE_MPI
|
||
pmesh = dynamic_cast<ParMesh *>(mesh);
|
||
if (pmesh) { parallel = true; }
|
||
#endif
|
||
}
|
||
void SetOrder(int order_) { order = order_ ; }
|
||
void SetDiffusionCoefficient(Coefficient * diffcf_) { diffcf = diffcf_; }
|
||
void SetMassCoefficient(Coefficient * masscf_) { masscf = masscf_; }
|
||
void SetRHSCoefficient(Coefficient * rhscf_) { rhscf = rhscf_; }
|
||
void SetEssentialBoundary(const Array<int> & ess_bdr_) { ess_bdr = ess_bdr_;};
|
||
void SetNeumannBoundary(const Array<int> & neumann_bdr_) { neumann_bdr = neumann_bdr_;};
|
||
void SetNeumannData(Coefficient * neumann_cf_) {neumann_cf = neumann_cf_;}
|
||
void SetEssBdrData(Coefficient * essbdr_cf_) {essbdr_cf = essbdr_cf_;}
|
||
void SetGradientData(VectorCoefficient * gradient_cf_) {gradient_cf = gradient_cf_;}
|
||
|
||
void ResetFEM();
|
||
void SetupFEM();
|
||
|
||
void UpdateEssentialTDofs();
|
||
void AssembleDiffusionBilinear(bool update_ess_tdofs=true);
|
||
void Solve();
|
||
GridFunction * GetFEMSolution();
|
||
LinearForm * GetLinearForm() {return b;}
|
||
#ifdef MFEM_USE_MPI
|
||
ParGridFunction * GetParFEMSolution();
|
||
ParLinearForm * GetParLinearForm()
|
||
{
|
||
if (parallel)
|
||
{
|
||
return dynamic_cast<ParLinearForm *>(b);
|
||
}
|
||
else
|
||
{
|
||
MFEM_ABORT("Wrong code path. Call GetLinearForm");
|
||
return nullptr;
|
||
}
|
||
}
|
||
#endif
|
||
|
||
~DiffusionSolver();
|
||
|
||
};
|
||
|
||
/**
|
||
* @brief Class for solving linear elasticity:
|
||
*
|
||
* -∇ ⋅ σ(u) = f in Ω + BCs
|
||
*
|
||
* where
|
||
*
|
||
* σ(u) = λ ∇⋅u I + μ (∇ u + ∇uᵀ)
|
||
*
|
||
*/
|
||
class LinearElasticitySolver
|
||
{
|
||
private:
|
||
Mesh * mesh = nullptr;
|
||
int order = 1;
|
||
Coefficient * lambda_cf = nullptr;
|
||
Coefficient * mu_cf = nullptr;
|
||
VectorCoefficient * essbdr_cf = nullptr;
|
||
VectorCoefficient * rhs_cf = nullptr;
|
||
|
||
// FEM solver
|
||
int dim;
|
||
FiniteElementCollection * fec = nullptr;
|
||
FiniteElementSpace * fes = nullptr;
|
||
Array<int> ess_bdr;
|
||
Array<int> neumann_bdr;
|
||
GridFunction * u = nullptr;
|
||
LinearForm * b = nullptr;
|
||
bool parallel;
|
||
#ifdef MFEM_USE_MPI
|
||
ParMesh * pmesh = nullptr;
|
||
ParFiniteElementSpace * pfes = nullptr;
|
||
#endif
|
||
|
||
public:
|
||
LinearElasticitySolver() { }
|
||
LinearElasticitySolver(Mesh * mesh_, int order_,
|
||
Coefficient * lambda_cf_, Coefficient * mu_cf_);
|
||
|
||
void SetMesh(Mesh * mesh_)
|
||
{
|
||
mesh = mesh_;
|
||
parallel = false;
|
||
#ifdef MFEM_USE_MPI
|
||
pmesh = dynamic_cast<ParMesh *>(mesh);
|
||
if (pmesh) { parallel = true; }
|
||
#endif
|
||
}
|
||
void SetOrder(int order_) { order = order_ ; }
|
||
void SetLameCoefficients(Coefficient * lambda_cf_, Coefficient * mu_cf_) { lambda_cf = lambda_cf_; mu_cf = mu_cf_; }
|
||
void SetRHSCoefficient(VectorCoefficient * rhs_cf_) { rhs_cf = rhs_cf_; }
|
||
void SetEssentialBoundary(const Array<int> & ess_bdr_) { ess_bdr = ess_bdr_;};
|
||
void SetNeumannBoundary(const Array<int> & neumann_bdr_) { neumann_bdr = neumann_bdr_;};
|
||
void SetEssBdrData(VectorCoefficient * essbdr_cf_) {essbdr_cf = essbdr_cf_;}
|
||
|
||
void ResetFEM();
|
||
void SetupFEM();
|
||
|
||
void Solve();
|
||
GridFunction * GetFEMSolution();
|
||
LinearForm * GetLinearForm() {return b;}
|
||
#ifdef MFEM_USE_MPI
|
||
ParGridFunction * GetParFEMSolution();
|
||
ParLinearForm * GetParLinearForm()
|
||
{
|
||
if (parallel)
|
||
{
|
||
return dynamic_cast<ParLinearForm *>(b);
|
||
}
|
||
else
|
||
{
|
||
MFEM_ABORT("Wrong code path. Call GetLinearForm");
|
||
return nullptr;
|
||
}
|
||
}
|
||
#endif
|
||
|
||
~LinearElasticitySolver();
|
||
|
||
};
|
||
|
||
/**
|
||
* @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
|
||
* ∫_Ω ρ dx = θ vol(Ω) as follows:
|
||
*
|
||
* 1. Compute the root of the R → R function
|
||
* f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
|
||
* using the Illinois method
|
||
* 2. Set ψ ← ψ + c.
|
||
*
|
||
* @param psi a GridFunction to be updated
|
||
* @param alpha_grad alpha multiplied by gradient
|
||
* @param target_volume θ vol(Ω)
|
||
* @param tol Illinois iteration tolerance
|
||
* @param max_its Illinois maximum iteration number
|
||
* @return real_t Final volume (∫_Ω sigmoid(ψ) dx)
|
||
*/
|
||
real_t proj(GridFunction &psi, GridFunction &alpha_grad, real_t target_volume,
|
||
real_t tol = 1e-12, int max_its = 100)
|
||
{
|
||
#ifdef MFEM_USE_MPI
|
||
FiniteElementSpace *fes = psi.FESpace();
|
||
ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
|
||
#endif
|
||
ConstantCoefficient zero_cf(0.0);
|
||
real_t a = -alpha_grad.ComputeMaxError(zero_cf);
|
||
real_t b = -a;
|
||
real_t y = 0.0;
|
||
|
||
MappedGridFunctionCoefficient sigmoid_psi(
|
||
&psi, [&y](const real_t x) { return sigmoid(x + y); });
|
||
std::unique_ptr<LinearForm> int_sigmoid_psi;
|
||
#ifdef MFEM_USE_MPI
|
||
ParGridFunction *par_psi = dynamic_cast<ParGridFunction *>(&psi);
|
||
if (par_psi)
|
||
{
|
||
int_sigmoid_psi.reset(new ParLinearForm(par_psi->ParFESpace()));
|
||
}
|
||
else
|
||
{
|
||
int_sigmoid_psi.reset(new LinearForm(psi.FESpace()));
|
||
}
|
||
#else
|
||
int_sigmoid_psi.reset(new LinearForm(psi.FESpace()));
|
||
#endif
|
||
int_sigmoid_psi->AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
|
||
|
||
y = a;
|
||
int_sigmoid_psi->Assemble();
|
||
real_t f_a = int_sigmoid_psi->Sum(); // f_a := f(a) + θ vol(Ω)
|
||
|
||
y = b;
|
||
int_sigmoid_psi->Assemble();
|
||
real_t f_b = int_sigmoid_psi->Sum(); // f_b := f(b) + θ vol(Ω)
|
||
#ifdef MFEM_USE_MPI
|
||
if (pfes)
|
||
{
|
||
MPI_Allreduce(MPI_IN_PLACE, &f_a, 1, MPITypeMap<real_t>::mpi_type,
|
||
MPI_SUM, MPI_COMM_WORLD);
|
||
MPI_Allreduce(MPI_IN_PLACE, &f_b, 1, MPITypeMap<real_t>::mpi_type,
|
||
MPI_SUM, MPI_COMM_WORLD);
|
||
}
|
||
#endif
|
||
f_a -= target_volume; // f_a := f(a)
|
||
f_b -= target_volume; // f_b := f(b)
|
||
real_t c = 0.0;
|
||
real_t f_c = 0.0;
|
||
int side = 0;
|
||
|
||
bool done = false;
|
||
for (int k=0; k < max_its; k++)
|
||
{
|
||
c = (f_a * b - f_b * a) / (f_a - f_b);
|
||
|
||
if (abs(b - a) < tol * abs(b + a)) { done = true; break; }
|
||
|
||
y = c;
|
||
int_sigmoid_psi->Assemble();
|
||
f_c = int_sigmoid_psi->Sum(); // f_c := f(c) + θ vol(Ω)
|
||
#ifdef MFEM_USE_MPI
|
||
if (pfes)
|
||
{
|
||
MPI_Allreduce(MPI_IN_PLACE, &f_c, 1, MPITypeMap<real_t>::mpi_type,
|
||
MPI_SUM, MPI_COMM_WORLD);
|
||
}
|
||
#endif
|
||
f_c -= target_volume; // f_c := f(c)
|
||
|
||
if (f_c * f_b > 0)
|
||
{
|
||
b = c;
|
||
f_b = f_c;
|
||
if (side == -1) { f_a /= 2.0; }
|
||
side = -1;
|
||
}
|
||
else if (f_c * f_a > 0)
|
||
{
|
||
a = c;
|
||
f_a = f_c;
|
||
if (side == 1) { f_b /= 2.0; }
|
||
side = 1;
|
||
}
|
||
else
|
||
{
|
||
done = true; break;
|
||
}
|
||
}
|
||
if (!done)
|
||
{
|
||
mfem_warning("Projection reached maximum iteration without converging. "
|
||
"Result may not be accurate.");
|
||
}
|
||
y = 0.0;
|
||
psi += c;
|
||
int_sigmoid_psi->Assemble();
|
||
real_t material_volume = int_sigmoid_psi->Sum();
|
||
#ifdef MFEM_USE_MPI
|
||
if (pfes)
|
||
{
|
||
MPI_Allreduce(MPI_IN_PLACE, &material_volume, 1,
|
||
MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
|
||
}
|
||
#endif
|
||
return material_volume;
|
||
}
|
||
|
||
// Poisson solver
|
||
|
||
DiffusionSolver::DiffusionSolver(Mesh * mesh_, int order_,
|
||
Coefficient * diffcf_, Coefficient * rhscf_)
|
||
: mesh(mesh_), order(order_), diffcf(diffcf_), rhscf(rhscf_)
|
||
{
|
||
|
||
#ifdef MFEM_USE_MPI
|
||
pmesh = dynamic_cast<ParMesh *>(mesh);
|
||
if (pmesh) { parallel = true; }
|
||
#endif
|
||
|
||
SetupFEM();
|
||
}
|
||
|
||
void DiffusionSolver::SetupFEM()
|
||
{
|
||
dim = mesh->Dimension();
|
||
fec = new H1_FECollection(order, dim);
|
||
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
pfes = new ParFiniteElementSpace(pmesh, fec);
|
||
u = new ParGridFunction(pfes);
|
||
b = new ParLinearForm(pfes);
|
||
}
|
||
else
|
||
{
|
||
fes = new FiniteElementSpace(mesh, fec);
|
||
u = new GridFunction(fes);
|
||
b = new LinearForm(fes);
|
||
}
|
||
#else
|
||
fes = new FiniteElementSpace(mesh, fec);
|
||
u = new GridFunction(fes);
|
||
b = new LinearForm(fes);
|
||
#endif
|
||
*u=0.0;
|
||
|
||
if (!ess_bdr.Size())
|
||
{
|
||
if (mesh->bdr_attributes.Size())
|
||
{
|
||
ess_bdr.SetSize(mesh->bdr_attributes.Max());
|
||
ess_bdr = 1;
|
||
}
|
||
}
|
||
}
|
||
|
||
void DiffusionSolver::UpdateEssentialTDofs()
|
||
{
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
pfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
|
||
}
|
||
else
|
||
{
|
||
fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
|
||
}
|
||
#else
|
||
fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
|
||
#endif
|
||
}
|
||
|
||
void DiffusionSolver::AssembleDiffusionBilinear(bool update_ess_tdofs)
|
||
{
|
||
if (update_ess_tdofs)
|
||
{
|
||
UpdateEssentialTDofs();
|
||
}
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
a = new ParBilinearForm(pfes);
|
||
}
|
||
else
|
||
{
|
||
a = new BilinearForm(fes);
|
||
}
|
||
#else
|
||
a = new BilinearForm(fes);
|
||
#endif
|
||
a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
|
||
if (masscf)
|
||
{
|
||
a->AddDomainIntegrator(new MassIntegrator(*masscf));
|
||
}
|
||
a->Assemble();
|
||
a->FormSystemMatrix(ess_tdof_list, A);
|
||
}
|
||
|
||
void DiffusionSolver::Solve()
|
||
{
|
||
Vector B, X;
|
||
|
||
if (b)
|
||
{
|
||
delete b;
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
b = new ParLinearForm(pfes);
|
||
}
|
||
else
|
||
{
|
||
b = new LinearForm(fes);
|
||
}
|
||
#else
|
||
b = new LinearForm(fes);
|
||
#endif
|
||
}
|
||
if (rhscf)
|
||
{
|
||
b->AddDomainIntegrator(new DomainLFIntegrator(*rhscf));
|
||
}
|
||
if (neumann_cf)
|
||
{
|
||
MFEM_VERIFY(neumann_bdr.Size(), "neumann_bdr attributes not provided");
|
||
b->AddBoundaryIntegrator(new BoundaryLFIntegrator(*neumann_cf),neumann_bdr);
|
||
}
|
||
else if (gradient_cf)
|
||
{
|
||
MFEM_VERIFY(neumann_bdr.Size(), "neumann_bdr attributes not provided");
|
||
b->AddBoundaryIntegrator(new BoundaryNormalLFIntegrator(*gradient_cf),
|
||
neumann_bdr);
|
||
}
|
||
|
||
b->Assemble();
|
||
|
||
*u=0.0;
|
||
if (essbdr_cf)
|
||
{
|
||
u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
|
||
}
|
||
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
X.SetSize(pfes->TrueVSize());
|
||
B.SetSize(pfes->TrueVSize());
|
||
dynamic_cast<ParGridFunction*>(u)->ParallelAssemble(X);
|
||
dynamic_cast<ParLinearForm*>(b)->ParallelAssemble(B);
|
||
dynamic_cast<ParBilinearForm*>(a)->ParallelEliminateTDofsInRHS(
|
||
ess_tdof_list, X, B);
|
||
}
|
||
else
|
||
{
|
||
X.NewDataAndSize(u->GetData(), u->Size());
|
||
B.NewDataAndSize(b->GetData(), b->Size());
|
||
a->EliminateVDofsInRHS(ess_tdof_list, X, B);
|
||
}
|
||
#else
|
||
X.NewDataAndSize(u->GetData(), u->Size());
|
||
B.NewDataAndSize(b->GetData(), b->Size());
|
||
a->EliminateVDofsInRHS(ess_tdof_list, X, B);
|
||
#endif
|
||
|
||
CGSolver * cg = nullptr;
|
||
Solver * M = nullptr;
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
M = new HypreBoomerAMG;
|
||
dynamic_cast<HypreBoomerAMG*>(M)->SetPrintLevel(0);
|
||
cg = new CGSolver(pmesh->GetComm());
|
||
}
|
||
else
|
||
{
|
||
M = new GSSmoother((SparseMatrix&)(*A));
|
||
cg = new CGSolver;
|
||
}
|
||
#else
|
||
M = new GSSmoother((SparseMatrix&)(*A));
|
||
cg = new CGSolver;
|
||
#endif
|
||
cg->SetRelTol(1e-12);
|
||
cg->SetMaxIter(10000);
|
||
cg->SetPrintLevel(0);
|
||
cg->SetPreconditioner(*M);
|
||
cg->SetOperator(*A);
|
||
cg->Mult(B, X);
|
||
delete M;
|
||
delete cg;
|
||
a->RecoverFEMSolution(X, *b, *u);
|
||
}
|
||
|
||
GridFunction * DiffusionSolver::GetFEMSolution()
|
||
{
|
||
return u;
|
||
}
|
||
|
||
#ifdef MFEM_USE_MPI
|
||
ParGridFunction * DiffusionSolver::GetParFEMSolution()
|
||
{
|
||
if (parallel)
|
||
{
|
||
return dynamic_cast<ParGridFunction*>(u);
|
||
}
|
||
else
|
||
{
|
||
MFEM_ABORT("Wrong code path. Call GetFEMSolution");
|
||
return nullptr;
|
||
}
|
||
}
|
||
#endif
|
||
|
||
DiffusionSolver::~DiffusionSolver()
|
||
{
|
||
delete u; u = nullptr;
|
||
delete fes; fes = nullptr;
|
||
#ifdef MFEM_USE_MPI
|
||
delete pfes; pfes=nullptr;
|
||
#endif
|
||
delete fec; fec = nullptr;
|
||
delete b;
|
||
A.Clear();
|
||
delete a;
|
||
}
|
||
|
||
|
||
// Elasticity solver
|
||
|
||
LinearElasticitySolver::LinearElasticitySolver(Mesh * mesh_, int order_,
|
||
Coefficient * lambda_cf_, Coefficient * mu_cf_)
|
||
: mesh(mesh_), order(order_), lambda_cf(lambda_cf_), mu_cf(mu_cf_)
|
||
{
|
||
#ifdef MFEM_USE_MPI
|
||
pmesh = dynamic_cast<ParMesh *>(mesh);
|
||
if (pmesh) { parallel = true; }
|
||
#endif
|
||
SetupFEM();
|
||
}
|
||
|
||
void LinearElasticitySolver::SetupFEM()
|
||
{
|
||
dim = mesh->Dimension();
|
||
fec = new H1_FECollection(order, dim,BasisType::Positive);
|
||
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
pfes = new ParFiniteElementSpace(pmesh, fec, dim);
|
||
u = new ParGridFunction(pfes);
|
||
b = new ParLinearForm(pfes);
|
||
}
|
||
else
|
||
{
|
||
fes = new FiniteElementSpace(mesh, fec,dim);
|
||
u = new GridFunction(fes);
|
||
b = new LinearForm(fes);
|
||
}
|
||
#else
|
||
fes = new FiniteElementSpace(mesh, fec, dim);
|
||
u = new GridFunction(fes);
|
||
b = new LinearForm(fes);
|
||
#endif
|
||
*u=0.0;
|
||
|
||
if (!ess_bdr.Size())
|
||
{
|
||
if (mesh->bdr_attributes.Size())
|
||
{
|
||
ess_bdr.SetSize(mesh->bdr_attributes.Max());
|
||
ess_bdr = 1;
|
||
}
|
||
}
|
||
}
|
||
|
||
void LinearElasticitySolver::Solve()
|
||
{
|
||
GridFunction * x = nullptr;
|
||
OperatorPtr A;
|
||
Vector B, X;
|
||
Array<int> ess_tdof_list;
|
||
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
x = new ParGridFunction(pfes);
|
||
pfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
|
||
}
|
||
else
|
||
{
|
||
x = new GridFunction(fes);
|
||
fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
|
||
}
|
||
#else
|
||
x = new GridFunction(fes);
|
||
fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
|
||
#endif
|
||
*u=0.0;
|
||
if (b)
|
||
{
|
||
delete b;
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
b = new ParLinearForm(pfes);
|
||
}
|
||
else
|
||
{
|
||
b = new LinearForm(fes);
|
||
}
|
||
#else
|
||
b = new LinearForm(fes);
|
||
#endif
|
||
}
|
||
if (rhs_cf)
|
||
{
|
||
b->AddDomainIntegrator(new VectorDomainLFIntegrator(*rhs_cf));
|
||
}
|
||
|
||
b->Assemble();
|
||
|
||
*x = 0.0;
|
||
|
||
BilinearForm * a = nullptr;
|
||
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
a = new ParBilinearForm(pfes);
|
||
}
|
||
else
|
||
{
|
||
a = new BilinearForm(fes);
|
||
}
|
||
#else
|
||
a = new BilinearForm(fes);
|
||
#endif
|
||
a->AddDomainIntegrator(new ElasticityIntegrator(*lambda_cf, *mu_cf));
|
||
a->Assemble();
|
||
if (essbdr_cf)
|
||
{
|
||
u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
|
||
}
|
||
a->FormLinearSystem(ess_tdof_list, *x, *b, A, X, B);
|
||
|
||
CGSolver * cg = nullptr;
|
||
Solver * M = nullptr;
|
||
#ifdef MFEM_USE_MPI
|
||
if (parallel)
|
||
{
|
||
M = new HypreBoomerAMG;
|
||
dynamic_cast<HypreBoomerAMG*>(M)->SetPrintLevel(0);
|
||
cg = new CGSolver(pmesh->GetComm());
|
||
}
|
||
else
|
||
{
|
||
M = new GSSmoother((SparseMatrix&)(*A));
|
||
cg = new CGSolver;
|
||
}
|
||
#else
|
||
M = new GSSmoother((SparseMatrix&)(*A));
|
||
cg = new CGSolver;
|
||
#endif
|
||
cg->SetRelTol(1e-10);
|
||
cg->SetMaxIter(10000);
|
||
cg->SetPrintLevel(0);
|
||
cg->SetPreconditioner(*M);
|
||
cg->SetOperator(*A);
|
||
cg->Mult(B, X);
|
||
delete M;
|
||
delete cg;
|
||
a->RecoverFEMSolution(X, *b, *x);
|
||
*u+=*x;
|
||
delete a;
|
||
delete x;
|
||
}
|
||
|
||
GridFunction * LinearElasticitySolver::GetFEMSolution()
|
||
{
|
||
return u;
|
||
}
|
||
|
||
#ifdef MFEM_USE_MPI
|
||
ParGridFunction * LinearElasticitySolver::GetParFEMSolution()
|
||
{
|
||
if (parallel)
|
||
{
|
||
return dynamic_cast<ParGridFunction*>(u);
|
||
}
|
||
else
|
||
{
|
||
MFEM_ABORT("Wrong code path. Call GetFEMSolution");
|
||
return nullptr;
|
||
}
|
||
}
|
||
#endif
|
||
|
||
LinearElasticitySolver::~LinearElasticitySolver()
|
||
{
|
||
delete u; u = nullptr;
|
||
delete fes; fes = nullptr;
|
||
#ifdef MFEM_USE_MPI
|
||
delete pfes; pfes=nullptr;
|
||
#endif
|
||
delete fec; fec = nullptr;
|
||
delete b;
|
||
}
|
||
|
||
} // namespace mfem
|
||
|