71 const double small = 0.001, big = 0.01;
76 const double X = x(0), Y = x(1);
77 const double ind = std::tanh((10*(Y-0.5) + std::sin(4.0*M_PI*X)) + 1) -
78 std::tanh((10*(Y-0.5) + std::sin(4.0*M_PI*X)) - 1);
80 return ind * small + (1.0 - ind) * big;
87 const double xc = x(0) - 0.5, yc = x(1) - 0.5;
88 const double r = sqrt(xc*xc + yc*yc);
89 double r1 = 0.15;
double r2 = 0.35;
double sf=30.0;
90 val = 0.5*(std::tanh(sf*(r-r1)) - std::tanh(sf*(r-r2)));
91 if (val > 1.) {val = 1;}
93 return val * small + (1.0 - val) * big;
99 const double X = x(0), Y = x(1);
100 const double r1 = 0.45, r2 = 0.55;
101 const double sf = 40.0;
103 double val = 0.5 * ( std::tanh(sf*(X-r1)) - std::tanh(sf*(X-r2)) +
104 std::tanh(sf*(Y-r1)) - std::tanh(sf*(Y-r2)) );
105 if (val > 1.) { val = 1.0; }
107 return val * small + (1.0 - val) * big;
113 double r1,r2,val,rval;
117 r1= 0.25; r2 = 0.25; rval = 0.1;
118 double xc = x(0) - r1, yc = x(1) - r2;
119 double r = sqrt(xc*xc+yc*yc);
120 val = 0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
124 xc = x(0) - r1, yc = x(1) - r2;
125 r = sqrt(xc*xc+yc*yc);
126 val += (0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
130 xc = x(0) - r1, yc = x(1) - r2;
131 r = sqrt(xc*xc+yc*yc);
132 val += 0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
136 xc = x(0) - r1, yc = x(1) - r2;
137 r = sqrt(xc*xc+yc*yc);
138 val += 0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*(r-rval)));
139 if (val > 1.0) {val = 1.;}
140 if (val < 0.0) {val = 0.;}
142 return val * small + (1.0 - val) * big;
149 double X = x(0)-0.5, Y = x(1)-0.5;
150 double rval = std::sqrt(X*X + Y*Y);
151 double thval = 60.*M_PI/180.;
153 Xmod = X*std::cos(thval) + Y*std::sin(thval);
154 Ymod= -X*std::sin(thval) + Y*std::cos(thval);
155 X = Xmod+0.5; Y = Ymod+0.5;
156 double r1 = 0.45;
double r2 = 0.55;
double sf=30.0;
157 val = ( 0.5*(1+std::tanh(sf*(X-r1))) - 0.5*(1+std::tanh(sf*(X-r2)))
158 + 0.5*(1+std::tanh(sf*(Y-r1))) - 0.5*(1+std::tanh(sf*(Y-r2))) );
159 if (rval > 0.4) {val = 0.;}
160 if (val > 1.0) {val = 1.;}
161 if (val < 0.0) {val = 0.;}
163 return val * small + (1.0 - val) * big;
169 const double xc = x(0) - 0.0, yc = x(1) - 0.5;
170 const double r = sqrt(xc*xc + yc*yc);
171 double r1 = 0.45;
double r2 = 0.55;
double sf=30.0;
172 val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
173 if (val > 1.) {val = 1;}
174 if (val < 0.) {val = 0;}
176 return val * small + (1.0 - val) * big;
188 HessianCoefficient(
int dim,
int type_)
199 K(0, 0) = 1.0 + 3.0 * std::sin(M_PI*pos(0));
206 const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
207 const double r = sqrt(xc*xc + yc*yc);
208 double r1 = 0.15;
double r2 = 0.35;
double sf=30.0;
209 const double eps = 0.5;
211 const double tan1 = std::tanh(sf*(r-r1)),
212 tan2 = std::tanh(sf*(r-r2));
214 K(0, 0) = eps + 1.0 * (tan1 - tan2);
226 int main (
int argc,
char *argv[])
230 MPI_Init(&argc, &argv);
231 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
232 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
235 const char *mesh_file =
"icf.mesh";
236 int mesh_poly_deg = 1;
242 double lim_const = 0.0;
245 int newton_iter = 10;
246 double newton_rtol = 1e-12;
248 int max_lin_iter = 100;
249 bool move_bnd =
true;
251 bool normalization =
false;
252 bool visualization =
true;
253 int verbosity_level = 0;
257 args.
AddOption(&mesh_file,
"-m",
"--mesh",
258 "Mesh file to use.");
259 args.
AddOption(&mesh_poly_deg,
"-o",
"--order",
260 "Polynomial degree of mesh finite element space.");
261 args.
AddOption(&rs_levels,
"-rs",
"--refine-serial",
262 "Number of times to refine the mesh uniformly in serial.");
263 args.
AddOption(&rp_levels,
"-rp",
"--refine-parallel",
264 "Number of times to refine the mesh uniformly in parallel.");
265 args.
AddOption(&jitter,
"-ji",
"--jitter",
266 "Random perturbation scaling factor.");
267 args.
AddOption(&metric_id,
"-mid",
"--metric-id",
268 "Mesh optimization metric:\n\t"
269 "1 : |T|^2 -- 2D shape\n\t"
270 "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
271 "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
272 "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
273 "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
274 "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
275 "55 : (tau-1)^2 -- 2D size\n\t"
276 "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
277 "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
278 "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
279 "211: (tau-1)^2-tau+sqrt(tau^2) -- 2D untangling\n\t"
280 "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
281 "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
282 "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
283 "303: (|T|^2)/3*tau^(2/3)-1 -- 3D shape\n\t"
284 "315: (tau-1)^2 -- 3D size\n\t"
285 "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D size\n\t"
286 "321: |T-T^-t|^2 -- 3D shape+size\n\t"
287 "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling");
288 args.
AddOption(&target_id,
"-tid",
"--target-id",
289 "Target (ideal element) type:\n\t"
290 "1: Ideal shape, unit size\n\t"
291 "2: Ideal shape, equal size\n\t"
292 "3: Ideal shape, initial size\n\t"
293 "4: Given full analytic Jacobian (in physical space)\n\t"
294 "5: Ideal shape, given size (in physical space)");
295 args.
AddOption(&lim_const,
"-lc",
"--limit-const",
"Limiting constant.");
296 args.
AddOption(&quad_type,
"-qt",
"--quad-type",
297 "Quadrature rule type:\n\t"
298 "1: Gauss-Lobatto\n\t"
299 "2: Gauss-Legendre\n\t"
300 "3: Closed uniform points");
301 args.
AddOption(&quad_order,
"-qo",
"--quad_order",
302 "Order of the quadrature rule.");
303 args.
AddOption(&newton_iter,
"-ni",
"--newton-iters",
304 "Maximum number of Newton iterations.");
305 args.
AddOption(&newton_rtol,
"-rtol",
"--newton-rel-tolerance",
306 "Relative tolerance for the Newton solver.");
307 args.
AddOption(&lin_solver,
"-ls",
"--lin-solver",
308 "Linear solver: 0 - l1-Jacobi, 1 - CG, 2 - MINRES.");
309 args.
AddOption(&max_lin_iter,
"-li",
"--lin-iter",
310 "Maximum number of iterations in the linear solve.");
311 args.
AddOption(&move_bnd,
"-bnd",
"--move-boundary",
"-fix-bnd",
313 "Enable motion along horizontal and vertical boundaries.");
314 args.
AddOption(&combomet,
"-cmb",
"--combo-met",
"-no-cmb",
"--no-combo-met",
315 "Combination of metrics.");
316 args.
AddOption(&normalization,
"-nor",
"--normalization",
"-no-nor",
317 "--no-normalization",
318 "Make all terms in the optimization functional unitless.");
319 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
320 "--no-visualization",
321 "Enable or disable GLVis visualization.");
322 args.
AddOption(&verbosity_level,
"-vl",
"--verbosity-level",
323 "Set the verbosity level - 0, 1, or 2.");
333 Mesh *mesh =
new Mesh(mesh_file, 1, 1,
false);
338 cout <<
"Mesh curvature: ";
340 else { cout <<
"(NONE)"; }
353 if (mesh_poly_deg <= 0)
382 double vol_loc = 0.0;
384 for (
int i = 0; i < pmesh->
GetNE(); i++)
390 for (
int j = 0; j < dofs.
Size(); j++)
392 h0(dofs[j]) = min(h0(dofs[j]), hi);
397 MPI_Allreduce(&vol_loc, &volume, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
398 const double small_phys_size = pow(volume, 1.0 / dim) / 100.0;
410 for (
int i = 0; i < pfespace->
GetNDofs(); i++)
412 for (
int d = 0; d <
dim; d++)
418 for (
int i = 0; i < pfespace->
GetNBE(); i++)
423 for (
int j = 0; j < vdofs.
Size(); j++) { rdm(vdofs[j]) = 0.0; }
434 ostringstream mesh_name;
435 mesh_name <<
"perturbed." << setfill(
'0') << setw(6) << myid;
436 ofstream mesh_ofs(mesh_name.str().c_str());
437 mesh_ofs.precision(8);
438 pmesh->
Print(mesh_ofs);
446 double tauval = -0.1;
470 if (myid == 0) { cout <<
"Unknown metric_id: " << metric_id << endl; }
475 HessianCoefficient *adapt_coeff = NULL;
488 adapt_coeff =
new HessianCoefficient(dim, 1);
505 if (myid == 0) { cout <<
"Unknown target_id: " << target_id << endl; }
509 if (target_c == NULL)
522 case 2: ir = &
IntRules.
Get(geom_type, quad_order);
break;
525 if (myid == 0) { cout <<
"Unknown quad_type: " << quad_type << endl; }
529 { cout <<
"Quadrature points per cell: " << ir->
GetNPoints() << endl; }
539 if (normalization) { dist = small_phys_size; }
541 if (lim_const != 0.0) { he_nlf_integ->
EnableLimiting(x0, dist, lim_coeff); }
559 if (normalization) { MFEM_ABORT(
"Not implemented."); }
587 char title[] =
"Initial metric values";
596 if (move_bnd ==
false)
606 for (
int i = 0; i < pmesh->
GetNBE(); i++)
609 MFEM_VERIFY(!(dim == 2 && attr == 3),
610 "Boundary attribute 3 must be used only for 3D meshes. "
611 "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
612 "components, rest for free nodes), or use -fix-bnd.");
613 if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
614 if (attr == 4) { n += nd *
dim; }
618 for (
int i = 0; i < pmesh->
GetNBE(); i++)
624 for (
int j = 0; j < nd; j++)
625 { ess_vdofs[n++] = vdofs[j]; }
629 for (
int j = 0; j < nd; j++)
630 { ess_vdofs[n++] = vdofs[j+nd]; }
634 for (
int j = 0; j < nd; j++)
635 { ess_vdofs[n++] = vdofs[j+2*nd]; }
639 for (
int j = 0; j < vdofs.
Size(); j++)
640 { ess_vdofs[n++] = vdofs[j]; }
649 const double linsol_rtol = 1e-12;
654 else if (lin_solver == 1)
675 const int NE = pmesh->
GetNE();
676 for (
int i = 0; i < NE; i++)
686 MPI_Allreduce(&tauval, &minJ0, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
689 { cout <<
"Minimum det(J) of the original mesh is " << tauval << endl; }
703 { cout <<
"TMOPNewtonSolver is used (as all det(J) > 0)." << endl; }
707 if ( (dim == 2 && metric_id != 22 && metric_id != 252) ||
708 (dim == 3 && metric_id != 352) )
711 { cout <<
"The mesh is inverted. Use an untangling metric.\n"; }
714 double h0min = h0.Min(), h0min_all;
715 MPI_Allreduce(&h0min, &h0min_all, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
716 tauval -= 0.01 * h0min_all;
719 { cout <<
"TMOPDescentNewtonSolver is used (as some det(J) < 0).\n"; }
731 cout <<
"NewtonIteration: rtol = " << newton_rtol <<
" not achieved."
739 ostringstream mesh_name;
740 mesh_name <<
"optimized." << setfill(
'0') << setw(6) << myid;
741 ofstream mesh_ofs(mesh_name.str().c_str());
742 mesh_ofs.precision(8);
743 pmesh->
Print(mesh_ofs);
748 double metric_part = fin_energy;
749 if (lim_const != 0.0)
757 cout <<
"Initial strain energy: " << init_energy
758 <<
" = metrics: " << init_energy
759 <<
" + limiting term: " << 0.0 << endl;
760 cout <<
" Final strain energy: " << fin_energy
761 <<
" = metrics: " << metric_part
762 <<
" + limiting term: " << fin_energy - metric_part << endl;
763 cout <<
"The strain energy decreased by: " << setprecision(12)
764 << (init_energy - fin_energy) * 100.0 / init_energy <<
" %." << endl;
770 char title[] =
"Final metric values";
781 sock.
open(
"localhost", 19916);
782 sock <<
"solution\n";
788 sock <<
"window_title 'Displacements'\n"
789 <<
"window_geometry "
790 << 1200 <<
" " << 0 <<
" " << 600 <<
" " << 600 <<
"\n"
791 <<
"keys jRmclA" << endl;
814 const double r = sqrt(x(0)*x(0) + x(1)*x(1) + 1e-12);
815 const double den = 0.002;
816 double l2 = 0.2 + 0.5 * (std::tanh((r-0.16)/den) - std::tanh((r-0.17)/den)
817 + std::tanh((r-0.23)/den) - std::tanh((r-0.24)/den));
int GetNPoints() const
Returns the number of the points in the integration rule.
void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
virtual void SetParDiscreteTargetSpec(ParGridFunction &tspec)
int Size() const
Logical size of the array.
Shifted barrier form of metric 56 (area, ideal barrier metric), 2D.
Shifted barrier form of 3D metric 16 (volume, ideal barrier metric), 3D.
Conjugate gradient method.
int GetNDofs() const
Returns number of degrees of freedom.
Class for an integration rule - an Array of IntegrationPoint.
int DofToVDof(int dof, int vd, int ndofs=-1) const
Shape & volume, ideal barrier metric, 3D.
Data type for scaled Jacobi-type smoother of sparse matrix.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Subclass constant coefficient.
int GetNBE() const
Returns number of boundary elements.
Shape, ideal barrier metric, 3D.
Allows negative Jacobians. Used for untangling.
void SaveAsOne(std::ostream &out=mfem::out)
Merge the local grid functions.
Container class for integration rules.
Data type dense matrix using column-major storage.
Shape, ideal barrier metric, 3D.
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Shape, ideal barrier metric, 2D.
Area, ideal barrier metric, 2D.
int GetNE() const
Returns number of elements.
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the ParGridFunction.
Abstract parallel finite element space.
virtual void ProjectCoefficient(Coefficient &coeff)
void Randomize(int seed=0)
Set random values in the vector.
int main(int argc, char *argv[])
double weight_fun(const Vector &x)
Shape & area, ideal barrier metric, 2D.
Shifted barrier form of metric 2 (shape, ideal barrier metric), 2D.
Shape, ideal barrier metric, 2D.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
void SetPrintLevel(int print_lvl)
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Adds a limiting term to the integrator (general version).
int GetNBE() const
Returns number of boundary elements in the mesh.
void SetCoefficient(Coefficient &w1)
Sets a scaling Coefficient for the quality metric term of the integrator.
Volume, ideal barrier metric, 3D.
void PrintAsOne(std::ostream &out=mfem::out)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto)
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
void SetMaxIter(int max_it)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Newton's method for solving F(x)=b for a given operator F.
virtual void Print(std::ostream &out=mfem::out) const
Version of QuadraticFECollection with positive basis functions.
int GetAttribute() const
Return element's attribute.
double ind_values(const Vector &x)
void SetNodalFESpace(FiniteElementSpace *nfes)
void PrintUsage(std::ostream &out) const
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
double GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
void SetAbsTol(double atol)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void SetRelTol(double rtol)
Shape, ideal barrier metric, 2D.
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Area, ideal barrier metric, 2D.
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
virtual void SetAnalyticTargetSpec(Coefficient *sspec, VectorCoefficient *vspec, MatrixCoefficient *mspec)
Shape & area metric, 2D.
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
Class for integration point with weight.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i'th element.
void PrintOptions(std::ostream &out) const
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
int open(const char hostname[], int port)
class for C-function coefficient
IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform)
void GetNodes(Vector &node_coord) const
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Arbitrary order H1-conforming (continuous) finite elements.
TargetType
Target-matrix construction algorithms supported by this class.
void SetDiscreteAdaptTC(DiscreteAdaptTC *tc)
Class for parallel grid function.
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement for the i'th boundary element.
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th boundary element.
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Class for parallel meshes.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
void ParEnableNormalization(const ParGridFunction &x)
double GetElementVolume(int i)
const Element * GetBdrElement(int i) const
Shape, ideal barrier metric, 3D.
Metric without a type, 2D.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.