33 MFEM_VERIFY(nodes0.
Size() == new_mesh_nodes.
Size(),
34 "AdvectorCG assumes fixed mesh topology!");
41 ncomp =
space->GetVDim();
43 const int dof_cnt = field0.
Size() / ncomp;
47 for (
int i = 0; i < ncomp; i++)
51 new_field_temp.
MakeRef(new_field, i*dof_cnt, dof_cnt);
55 new_field_temp.
SetSize(dof_cnt);
56 for (
int j = 0; j < dof_cnt; j++)
58 new_field_temp(j) = new_field(i + j*ncomp);
61 ComputeAtNewPositionScalar(new_mesh_nodes, new_field_temp);
64 for (
int j = 0; j < dof_cnt; j++)
66 new_field(i + j*ncomp) = new_field_temp(j);
72 nodes0 = new_mesh_nodes;
75void AdvectorCG::ComputeAtNewPositionScalar(
const Vector &new_mesh_nodes,
83 MFEM_VERIFY(m != NULL,
"No mesh has been given to the AdaptivityEvaluator.");
112 MFEM_VERIFY(oper != NULL,
113 "No FE space has been given to the AdaptivityEvaluator.");
114 ode_solver.
Init(*oper);
117 real_t h_min = std::numeric_limits<real_t>::infinity();
118 for (
int i = 0; i < m->
GetNE(); i++)
126 for (
int i = 0; i < s; i++)
131 vel +=
u(i+j*s)*
u(i+j*s);
133 v_max = std::max(v_max,
vel);
139 real_t v_loc = v_max, h_loc = h_min;
157 v_max = std::sqrt(v_max);
158 real_t dt = dt_scale * h_min / v_max;
161 bool last_step =
false;
169 ode_solver.
Step(new_field, t, dt);
186 for (
int i = 0; i < new_field.
Size(); i++)
188 if (new_field(i) < glob_minv) { new_field(i) = glob_minv; }
189 if (new_field(i) > glob_maxv) { new_field(i) = glob_maxv; }
204 x0(x_start), x_now(*fes.
GetMesh()->GetNodes()),
205 u(
vel), u_coeff(&
u), M(&fes), K(&fes), al(al)
228 K.BilinearForm::operator=(0.0);
232 M.BilinearForm::operator=(0.0);
235#ifdef MFEM_USE_SINGLE
238 const real_t rtol = 1e-12;
251 PCG(*A, S, B, X, 0, 100, rtol, 0.0);
258 PCG(
M.
SpMat(), S, rhs, di_dt, 0, 100, rtol, 0.0);
269 x0(x_start), x_now(*pfes.
GetMesh()->GetNodes()),
270 u(
vel), u_coeff(&
u), M(&pfes), K(&pfes), al(al)
301 K.BilinearForm::operator=(0.0);
305 M.BilinearForm::operator=(0.0);
330#ifdef MFEM_USE_SINGLE
338 lin_solver.
Mult(*RHS, X);
359 if (m->
GetNodes()->FESpace()->IsDGSpace())
361 MFEM_ABORT(
"InterpolatorFP is not supported for periodic meshes yet.");
364 const real_t rel_bbox_el = 0.1;
365 const real_t newton_tol = 1.0e-12;
366 const int npts_at_once = 256;
380 finder->
Setup(*m, rel_bbox_el, newton_tol, npts_at_once);
383 field0_gf = init_field;
387 Vector &new_field,
int nodes_ordering)
392 if (fes_new_field ==
nullptr && new_mesh_nodes.
Size() != nodes0.
Size())
394 MFEM_WARNING(
"Deprecated -- use ComputeAtGivenPositions() instead!");
400 (fes_new_field) ? fes_new_field : field0_gf.
FESpace();
408 finder->
Interpolate(mapped_nodes, field0_gf, new_field);
412 finder->
Interpolate(new_mesh_nodes, field0_gf, new_field, nodes_ordering);
417 Vector &values,
int p_ordering)
419 finder->
Interpolate(positions, field0_gf, values, p_ordering);
431 MFEM_VERIFY(!(
parallel && p_nlf == NULL),
"Invalid Operator subclass.");
440 MFEM_VERIFY(!(serial && nlf == NULL),
"Invalid Operator subclass.");
453 if (!cP) { d_loc = d_in; }
454 else { cP->
Mult(d_in, d_loc); }
465 real_t init_fit_avg_err, init_fit_max_err = 0.0;
474 mfem::out <<
"TMOPNewtonSolver converged "
475 "based on the surface fitting error.\n";
486 mfem::out <<
"TMOPNewtonSolver terminated "
487 "based on max number of times surface fitting weight can"
497 const bool untangling = (min_detT_in <= 0.0) ?
true :
false;
498 const real_t untangle_factor = 1.5;
503 MFEM_VERIFY(
min_det_ptr != NULL,
" Initial mesh was valid, but"
504 " intermediate mesh is invalid. Contact TMOP Developers.");
506 "This setup is not supported. Contact TMOP Developers.");
510 const bool have_b = (
b.Size() ==
Height());
513 bool x_out_ok =
false;
514 real_t energy_out = 0.0, min_detT_out;
516 real_t avg_fit_err, max_fit_err = 0.0;
525 for (
int i = 0; i < 12; i++)
535 add(d_in, -scale,
c, d_out);
539 if (!cP) { d_loc = d_out; }
540 else { cP->
Mult(d_out, d_loc); }
553 mfem::out <<
"Scale = " << scale <<
" Neg det(J) found.\n";
555 scale *= detJ_factor;
continue;
557 if (untangling ==
true && min_detT_out < *
min_det_ptr)
562 mfem::out <<
"Scale = " << scale <<
" Neg det(J) decreased.\n";
564 scale *= detJ_factor;
continue;
570 if (untangling) { x_out_ok =
true;
break; }
580 if (max_fit_err >= 1.2*init_fit_max_err)
584 mfem::out <<
"Scale = " << scale <<
" Surf fit err increased.\n";
586 scale *= 0.5;
continue;
601 if (energy_out > energy_in + 0.2*fabs(energy_in) ||
602 std::isnan(energy_out) != 0)
606 mfem::out <<
"Scale = " << scale <<
" Increasing energy: "
607 << energy_in <<
" --> " << energy_out <<
'\n';
609 scale *= 0.5;
continue;
614 if (have_b) {
r -=
b; }
617 if (norm_out > 1.2*norm_in)
621 mfem::out <<
"Scale = " << scale <<
" Norm increased: "
622 << norm_in <<
" --> " << norm_out <<
'\n';
624 scale *= 0.5;
continue;
626 else { x_out_ok =
true;
break; }
632 if (min_detT_out > 0.0)
637 {
mfem::out <<
"The mesh has been untangled at the used points!\n"; }
639 else { *
min_det_ptr = untangle_factor * min_detT_out; }
648 << min_detT_in <<
" -> " << min_detT_out
649 <<
" with " << scale <<
" scaling.\n";
654 << energy_in <<
" --> " << energy_out <<
" or "
655 << (energy_in - energy_out) / energy_in * 100.0
656 <<
"% with " << scale <<
" scaling.\n";
660 if (x_out_ok ==
false) { scale = 0.0; }
673 const Operator *P = fes_mesh_nodes->GetProlongationMatrix();
675 periodic = fes_mesh_nodes->IsDGSpace();
679 "The input's size must be the tdof size of the mesh nodes.");
685 "The input's size must match the size of the mesh nodes.");
691 for (
int i = 0; i < integs.
Size(); i++)
703 else { MFEM_ABORT(
"Invalid solver_type"); }
710 if (Pd) { Pd->
Mult(dx, dx_loc); }
711 else { dx_loc = dx; }
718 for (
int i = 0; i < integs.
Size(); i++)
733 for (
int i = 0; i < integs.
Size(); i++)
744 for (
int j = 0; j < ati.
Size(); j++)
746 ati[j]->UpdateSurfaceFittingWeight(factor);
761 for (
int i = 0; i < integs.
Size(); i++)
773 for (
int j = 0; j < ati.
Size(); j++)
777 weight = ati[j]->GetSurfaceFittingWeight();
796 real_t err_avg_loc, err_max_loc;
797 for (
int i = 0; i < integs.
Size(); i++)
805 err_avg = std::max(err_avg_loc, err_avg);
806 err_max = std::max(err_max_loc, err_max);
813 for (
int j = 0; j < ati.
Size(); j++)
817 ati[j]->GetSurfaceFittingErrors(d_loc, err_avg_loc, err_max_loc);
818 err_avg = std::max(err_avg_loc, err_avg);
819 err_max = std::max(err_max_loc, err_max);
833 for (
int i = 0; i < integs.
Size(); i++)
847 for (
int j = 0; j < ati.
Size(); j++)
869 for (
int i = 0; i < integs.
Size(); i++)
881 for (
int j = 0; j < ati.
Size(); j++)
883 dtc = ati[j]->GetDiscreteAdaptTC();
896 else { dx_loc = dx; }
899 for (
int i = 0; i < integs.
Size(); i++)
914 for (
int j = 0; j < ati.
Size(); j++)
916 ati[j]->UpdateAfterMeshPositionChange(dx_loc, *dx_fes);
919 ati[j]->ComputeUntangleMetricQuantiles(dx_loc, *dx_fes);
939 mfem::out <<
"Avg/Max surface fitting error: " <<
942 mfem::out <<
"Min/Max surface fitting weight: " <<
943 fitweights.
Min() <<
" " << fitweights.
Max() <<
"\n";
980 if (
dim == 1 || mixed_mesh ||
983 for (
int i = 0; i < NE; i++)
1003 for (
int j = 0; j < nsp; j++)
1007 min_detJ = std::min(min_detJ, Jpr.
Det());
1021 MPI_MIN, p_nlf->ParFESpace()->GetComm());
1026 min_detJ /= Wideal.
Det();
1036 char *title,
int position)
1045 sock.
open(
"localhost", 19916);
1046 sock <<
"solution\n";
1052 sock <<
"window_title '"<< title <<
"'\n"
1053 <<
"window_geometry "
1054 << position <<
" " << 0 <<
" " << 600 <<
" " << 600 <<
"\n"
1064 char *title,
int position)
1071 sock <<
"solution\n";
1075 sock <<
"window_title '"<< title <<
"'\n"
1076 <<
"window_geometry "
1077 << position <<
" " << 0 <<
" " << 600 <<
" " << 600 <<
"\n"
1090 R_H1->
Mult(dx, dx_r);
1091 R_L2->AddMultTranspose(dx_r, x);
ParFiniteElementSpace * pfes
void SetInitialField(const Vector &init_nodes, const Vector &init_field) override
void ComputeAtNewPosition(const Vector &new_mesh_nodes, Vector &new_field, int nodes_ordering=Ordering::byNODES) override
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.
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
int Size() const
Return the logical size of the array.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
@ GaussLobatto
Closed type.
Conjugate gradient method.
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the Conjugate Gradient method.
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
Data type for scaled Jacobi-type smoother of sparse matrix.
Data type dense matrix using column-major storage.
real_t * Data() const
Returns the matrix data array.
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
void ResetUpdateFlags()
Used in combination with the Update methods to avoid extra computations.
FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary collection of points.
void Setup(Mesh &m, const double bb_t=0.1, const double newt_tol=1.0e-12, const int npt_max=256)
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
virtual const Operator * GetProlongationMatrix() const
void GetNodePositions(const Vector &mesh_nodes, Vector &fes_node_pos, int fes_nodes_ordering=Ordering::byNODES) const
Compute the space's node positions w.r.t. given mesh positions. The function uses FiniteElement::GetN...
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Ordering::Type GetOrdering() const
Return the ordering method.
int GetNE() const
Returns number of elements in the mesh.
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
const SparseMatrix * GetConformingProlongation() const
const FiniteElementCollection * FEColl() const
Mesh * GetMesh() const
Returns the mesh.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Class for grid function - Vector with associated FE space.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual void GetElementDofValues(int el, Vector &dof_vals) const
FiniteElementSpace * FESpace()
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Wrapper for hypre's parallel vector class.
Parallel smoothers in hypre.
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.
void ComputeAtGivenPositions(const Vector &positions, Vector &values, int p_ordering=Ordering::byNODES) override
Direct interpolation of field0_gf at the given positions.
void SetInitialField(const Vector &init_nodes, const Vector &init_field) override
void ComputeAtNewPosition(const Vector &new_mesh_nodes, Vector &new_field, int nodes_ordering=Ordering::byNODES) override
PrintLevel print_options
Output behavior for the iterative solver.
void SetRelTol(real_t rtol)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void SetMaxIter(int max_it)
void SetAbsTol(real_t atol)
real_t Norm(const Vector &x) const
Return the inner product norm of x, using the inner product defined by Dot()
Arbitrary order "L2-conforming" discontinuous finite elements.
void Mult(const Vector &b, Vector &x) const override
Solve the nonlinear system with right-hand side b.
void NodesUpdated()
This function should be called after the mesh node coordinates have been updated externally,...
Geometry::Type GetTypicalElementGeometry() const
If the local mesh is not empty, return GetElementGeometry(0); otherwise, return a typical Geometry pr...
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
int GetNE() const
Returns number of elements.
int Dimension() const
Dimension of the reference space used within the elements.
real_t GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
void GetNodes(Vector &node_coord) const
void SetNodes(const Vector &node_coord)
Updates the vertex/node locations. Invokes NodesUpdated().
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
void Mult(const Vector &b, Vector &x) const override
Solve the nonlinear system with right-hand side b.
Class for standard nodal finite elements.
void ReorderLexToNative(int ncomp, Vector &dofs) const
Pointer to an Operator of a specified type.
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Jacobi smoothing for a given bilinear form (no matrix necessary).
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int height
Dimension of the output / number of rows in the matrix.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Performs a single remap advection step in parallel.
ParAdvectorCGOper(const Vector &x_start, GridFunction &vel, ParFiniteElementSpace &pfes, AssemblyLevel al=AssemblyLevel::LEGACY, MemoryType mt=MemoryType::DEFAULT)
VectorGridFunctionCoefficient u_coeff
void Mult(const Vector &ind, Vector &di_dt) const override
Operator application: y=A(x).
Abstract parallel finite element space.
const Operator * GetProlongationMatrix() const override
ParMesh * GetParMesh() const
Class for parallel grid function.
void ParallelAssemble(Vector &tv) const
Returns the vector assembled on the true dofs.
void SaveAsOne(const char *fname, int precision=16) const
Class for parallel meshes.
void PrintAsOne(std::ostream &out=mfem::out, const std::string &comments="") const
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
void Step(Vector &x, real_t &t, real_t &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in].
Performs a single remap advection step in serial.
VectorGridFunctionCoefficient u_coeff
SerialAdvectorCGOper(const Vector &x_start, GridFunction &vel, FiniteElementSpace &fes, AssemblyLevel al=AssemblyLevel::LEGACY)
void Mult(const Vector &ind, Vector &di_dt) const override
Operator application: y=A(x).
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication.
void SetInitialMeshPos(const GridFunction *x0)
const Array< TMOP_Integrator * > & GetTMOPIntegrators() const
real_t surf_fit_max_err_limit
real_t ComputeScalingFactor(const Vector &d, const Vector &b) const override
virtual void GetSurfaceFittingError(const Vector &d_loc, real_t &err_avg, real_t &err_max) const
void Mult(const Vector &b, Vector &x) const override
Optimizes the mesh positions given by x.
real_t surf_fit_scale_factor
bool IsSurfaceFittingEnabled() const
Check if surface fitting is enabled.
real_t surf_fit_err_rel_change_limit
real_t MinDetJpr_3D(const FiniteElementSpace *, const Vector &) const
real_t surf_fit_weight_limit
int surf_fit_adapt_count_limit
const IntegrationRule & GetIntegrationRule(const FiniteElement &el) const
bool surf_fit_coeff_update
void UpdateSurfaceFittingWeight(real_t factor) const
Update surface fitting weight as surf_fit_weight *= factor.
real_t surf_fit_avg_err_prvs
real_t MinDetJpr_2D(const FiniteElementSpace *, const Vector &) const
void ProcessNewState(const Vector &dx) const override
bool compute_metric_quantile_flag
real_t ComputeMinDet(const Vector &d_loc, const FiniteElementSpace &fes) const
bool surf_fit_converge_error
void GetSurfaceFittingWeight(Array< real_t > &weights) const
Get the surface fitting weight for all the TMOP integrators.
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
void UpdateAfterMeshPositionChange(const Vector &d, const FiniteElementSpace &d_fes)
void GetSurfaceFittingErrors(const Vector &d_loc, real_t &err_avg, real_t &err_max)
void ComputeUntangleMetricQuantiles(const Vector &d, const FiniteElementSpace &fes)
real_t GetSurfaceFittingWeight()
Get the surface fitting weight.
void SetInitialMeshPos(const GridFunction *x0)
bool IsSurfaceFittingEnabled()
DiscreteAdaptTC * GetDiscreteAdaptTC() const
void UpdateSurfaceFittingWeight(real_t factor)
Update the surface fitting weight as surf_fit_coeff *= factor;.
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P....
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Base abstract class for first order time dependent operators.
virtual real_t GetTime() const
Read the currently set time.
real_t Max() const
Returns the maximal element of the vector.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
real_t Min() const
Returns the minimal element of the vector.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
real_t u(const Vector &xvec)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void add(const Vector &v1, const Vector &v2, Vector &v)
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, Mesh &mesh, char *title, int position)
void InterpolateTMOP_QualityMetric(TMOP_QualityMetric &metric, const TargetConstructor &tc, const Mesh &mesh, GridFunction &metric_gf)
Interpolates the metric's values at the nodes of metric_gf.
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
void GetPeriodicPositions(const Vector &x_0, const Vector &dx, const FiniteElementSpace &fesL2, const FiniteElementSpace &fesH1, Vector &x)
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
void subtract(const Vector &x, const Vector &y, Vector &z)
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
MemoryType
Memory types supported by MFEM.
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
std::function< real_t(const Vector &)> f(real_t mass_coeff)
void vel(const Vector &x, real_t t, Vector &u)
bool iterations
Detailed information about each iteration will be reported to mfem::out.
bool warnings
If a non-fatal problem has been detected some context-specific information will be reported to mfem::...
bool first_and_last
Information about the first and last iteration will be printed to mfem::out.
bool summary
A summary of the solver process will be reported after the last iteration to mfem::out.
Helper struct to convert a C++ type to an MPI type.