111 #include "../common/mfem-common.hpp" 116 using namespace mfem;
119 int main(
int argc,
char *argv[])
122 const char *mesh_file =
"icf.mesh";
123 int mesh_poly_deg = 1;
128 double lim_const = 0.0;
129 double adapt_lim_const = 0.0;
130 double surface_fit_const = 0.0;
134 int solver_iter = 20;
135 double solver_rtol = 1e-10;
136 int solver_art_type = 0;
138 int max_lin_iter = 100;
139 bool move_bnd =
true;
141 bool bal_expl_combo =
false;
142 bool hradaptivity =
false;
143 int h_metric_id = -1;
144 bool normalization =
false;
145 bool visualization =
true;
146 int verbosity_level = 0;
147 bool fdscheme =
false;
149 bool exactaction =
false;
150 const char *devopt =
"cpu";
154 bool surface_fit_adapt =
false;
155 double surface_fit_threshold = -10;
156 int mesh_node_ordering = 0;
157 int barrier_type = 0;
158 int worst_case_type = 0;
162 args.
AddOption(&mesh_file,
"-m",
"--mesh",
163 "Mesh file to use.");
164 args.
AddOption(&mesh_poly_deg,
"-o",
"--order",
165 "Polynomial degree of mesh finite element space.");
166 args.
AddOption(&rs_levels,
"-rs",
"--refine-serial",
167 "Number of times to refine the mesh uniformly in serial.");
168 args.
AddOption(&jitter,
"-ji",
"--jitter",
169 "Random perturbation scaling factor.");
170 args.
AddOption(&metric_id,
"-mid",
"--metric-id",
171 "Mesh optimization metric:\n\t" 173 "1 : |T|^2 -- 2D no type\n\t" 174 "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t" 175 "7 : |T-T^-t|^2 -- 2D shape+size\n\t" 176 "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t" 177 "14 : |T-I|^2 -- 2D shape+size+orientation\n\t" 178 "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t" 179 "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t" 180 "55 : (tau-1)^2 -- 2D size\n\t" 181 "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t" 182 "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t" 183 "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t" 184 "80 : (1-gamma)mu_2 + gamma mu_77 -- 2D shape+size\n\t" 185 "85 : |T-|T|/sqrt(2)I|^2 -- 2D shape+orientation\n\t" 186 "98 : (1/tau)|T-I|^2 -- 2D shape+size+orientation\n\t" 189 "301: (|T||T^-1|)/3-1 -- 3D shape\n\t" 190 "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t" 191 "303: (|T|^2)/3/tau^(2/3)-1 -- 3D shape\n\t" 192 "304: (|T|^3)/3^{3/2}/tau-1 -- 3D shape\n\t" 194 "313: (|T|^2)(tau-tau0)^(-2/3)/3 -- 3D untangling\n\t" 195 "315: (tau-1)^2 -- 3D no type\n\t" 196 "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D no type\n\t" 197 "321: |T-T^-t|^2 -- 3D shape+size\n\t" 198 "322: |T-adjT^-t|^2 -- 3D shape+size\n\t" 199 "323: |J|^3-3sqrt(3)ln(det(J))-3sqrt(3) -- 3D shape+size\n\t" 200 "328: (1-gamma) mu_301 + gamma mu_316 -- 3D shape+size\n\t" 201 "332: (1-gamma) mu_302 + gamma mu_315 -- 3D shape+size\n\t" 202 "333: (1-gamma) mu_302 + gamma mu_316 -- 3D shape+size\n\t" 203 "334: (1-gamma) mu_303 + gamma mu_316 -- 3D shape+size\n\t" 204 "347: (1-gamma) mu_304 + gamma mu_316 -- 3D shape+size\n\t" 206 "360: (|T|^3)/3^{3/2}-tau -- 3D shape\n\t" 208 "11 : (1/4*alpha)|A-(adjA)^T(W^TW)/omega|^2 -- 2D shape\n\t" 209 "36 : (1/alpha)|A-W|^2 -- 2D shape+size+orientation\n\t" 210 "107: (1/2*alpha)|A-|A|/|W|W|^2 -- 2D shape+orientation\n\t" 211 "126: (1-gamma)nu_11 + gamma*nu_14a -- 2D shape+size\n\t" 213 args.
AddOption(&target_id,
"-tid",
"--target-id",
214 "Target (ideal element) type:\n\t" 215 "1: Ideal shape, unit size\n\t" 216 "2: Ideal shape, equal size\n\t" 217 "3: Ideal shape, initial size\n\t" 218 "4: Given full analytic Jacobian (in physical space)\n\t" 219 "5: Ideal shape, given size (in physical space)");
220 args.
AddOption(&lim_const,
"-lc",
"--limit-const",
"Limiting constant.");
221 args.
AddOption(&adapt_lim_const,
"-alc",
"--adapt-limit-const",
222 "Adaptive limiting coefficient constant.");
223 args.
AddOption(&surface_fit_const,
"-sfc",
"--surface-fit-const",
224 "Surface preservation constant.");
225 args.
AddOption(&quad_type,
"-qt",
"--quad-type",
226 "Quadrature rule type:\n\t" 227 "1: Gauss-Lobatto\n\t" 228 "2: Gauss-Legendre\n\t" 229 "3: Closed uniform points");
230 args.
AddOption(&quad_order,
"-qo",
"--quad_order",
231 "Order of the quadrature rule.");
232 args.
AddOption(&solver_type,
"-st",
"--solver-type",
233 " Type of solver: (default) 0: Newton, 1: LBFGS");
234 args.
AddOption(&solver_iter,
"-ni",
"--newton-iters",
235 "Maximum number of Newton iterations.");
236 args.
AddOption(&solver_rtol,
"-rtol",
"--newton-rel-tolerance",
237 "Relative tolerance for the Newton solver.");
238 args.
AddOption(&solver_art_type,
"-art",
"--adaptive-rel-tol",
239 "Type of adaptive relative linear solver tolerance:\n\t" 240 "0: None (default)\n\t" 241 "1: Eisenstat-Walker type 1\n\t" 242 "2: Eisenstat-Walker type 2");
243 args.
AddOption(&lin_solver,
"-ls",
"--lin-solver",
248 "3: MINRES + Jacobi preconditioner\n\t" 249 "4: MINRES + l1-Jacobi preconditioner");
250 args.
AddOption(&max_lin_iter,
"-li",
"--lin-iter",
251 "Maximum number of iterations in the linear solve.");
252 args.
AddOption(&move_bnd,
"-bnd",
"--move-boundary",
"-fix-bnd",
254 "Enable motion along horizontal and vertical boundaries.");
255 args.
AddOption(&combomet,
"-cmb",
"--combo-type",
256 "Combination of metrics options:\n\t" 257 "0: Use single metric\n\t" 258 "1: Shape + space-dependent size given analytically\n\t" 259 "2: Shape + adapted size given discretely; shared target");
260 args.
AddOption(&bal_expl_combo,
"-bec",
"--balance-explicit-combo",
261 "-no-bec",
"--balance-explicit-combo",
262 "Automatic balancing of explicit combo metrics.");
263 args.
AddOption(&hradaptivity,
"-hr",
"--hr-adaptivity",
"-no-hr",
264 "--no-hr-adaptivity",
265 "Enable hr-adaptivity.");
266 args.
AddOption(&h_metric_id,
"-hmid",
"--h-metric",
267 "Same options as metric_id. Used to determine refinement" 268 " type for each element if h-adaptivity is enabled.");
269 args.
AddOption(&normalization,
"-nor",
"--normalization",
"-no-nor",
270 "--no-normalization",
271 "Make all terms in the optimization functional unitless.");
272 args.
AddOption(&fdscheme,
"-fd",
"--fd_approximation",
273 "-no-fd",
"--no-fd-approx",
274 "Enable finite difference based derivative computations.");
275 args.
AddOption(&exactaction,
"-ex",
"--exact_action",
276 "-no-ex",
"--no-exact-action",
277 "Enable exact action of TMOP_Integrator.");
278 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
279 "--no-visualization",
280 "Enable or disable GLVis visualization.");
281 args.
AddOption(&verbosity_level,
"-vl",
"--verbosity-level",
282 "Set the verbosity level - 0, 1, or 2.");
283 args.
AddOption(&adapt_eval,
"-ae",
"--adaptivity-evaluator",
284 "0 - Advection based (DEFAULT), 1 - GSLIB.");
285 args.
AddOption(&devopt,
"-d",
"--device",
286 "Device configuration string, see Device::Configure().");
287 args.
AddOption(&pa,
"-pa",
"--partial-assembly",
"-no-pa",
288 "--no-partial-assembly",
"Enable Partial Assembly.");
289 args.
AddOption(&n_hr_iter,
"-nhr",
"--n_hr_iter",
290 "Number of hr-adaptivity iterations.");
291 args.
AddOption(&n_h_iter,
"-nh",
"--n_h_iter",
292 "Number of h-adaptivity iterations per r-adaptivity" 294 args.
AddOption(&surface_fit_adapt,
"-sfa",
"--adaptive-surface-fit",
"-no-sfa",
295 "--no-adaptive-surface-fit",
296 "Enable or disable adaptive surface fitting.");
297 args.
AddOption(&surface_fit_threshold,
"-sft",
"--surf-fit-threshold",
298 "Set threshold for surface fitting. TMOP solver will" 299 "terminate when max surface fitting error is below this limit");
300 args.
AddOption(&mesh_node_ordering,
"-mno",
"--mesh_node_ordering",
301 "Ordering of mesh nodes." 302 "0 (default): byNodes, 1: byVDIM");
303 args.
AddOption(&barrier_type,
"-btype",
"--barrier-type",
305 "1 - Shifted Barrier," 306 "2 - Pseudo Barrier.");
307 args.
AddOption(&worst_case_type,
"-wctype",
"--worst-case-type",
319 if (h_metric_id < 0) { h_metric_id = metric_id; }
323 MFEM_VERIFY(strcmp(devopt,
"cpu")==0,
"HR-adaptivity is currently only" 324 " supported on cpus.");
330 Mesh *mesh =
new Mesh(mesh_file, 1, 1,
false);
341 if (mesh_poly_deg <= 0)
371 double mesh_volume = 0.0;
373 for (
int i = 0; i < mesh->
GetNE(); i++)
379 for (
int j = 0; j < dofs.
Size(); j++)
381 h0(dofs[j]) = min(h0(dofs[j]), hi);
385 const double small_phys_size = pow(mesh_volume, 1.0 /
dim) / 100.0;
398 for (
int i = 0; i < fespace->
GetNDofs(); i++)
400 for (
int d = 0; d <
dim; d++)
406 for (
int i = 0; i < fespace->
GetNBE(); i++)
411 for (
int j = 0; j < vdofs.
Size(); j++) { rdm(vdofs[j]) = 0.0; }
420 ofstream mesh_ofs(
"perturbed.mesh");
421 mesh->
Print(mesh_ofs);
429 double min_detJ = -0.1;
476 cout <<
"Unknown metric_id: " << metric_id << endl;
495 default: cout <<
"Metric_id not supported for h-adaptivity: " << h_metric_id <<
502 switch (barrier_type)
510 default: cout <<
"barrier_type not supported: " << barrier_type << endl;
515 switch (worst_case_type)
523 default: cout <<
"worst_case_type not supported: " << worst_case_type << endl;
528 if (barrier_type > 0 || worst_case_type > 0)
530 if (barrier_type > 0)
532 MFEM_VERIFY(metric_id == 4 || metric_id == 14 || metric_id == 66,
533 "Metric not supported for shifted/pseudo barriers.");
544 if (metric_id < 300 || h_metric_id < 300)
546 MFEM_VERIFY(
dim == 2,
"Incompatible metric for 3D meshes");
548 if (metric_id >= 300 || h_metric_id >= 300)
550 MFEM_VERIFY(
dim == 3,
"Incompatible metric for 2D meshes");
560 GridFunction size(&ind_fes), aspr(&ind_fes), ori(&ind_fes);
590 #ifdef MFEM_USE_GSLIB 593 MFEM_ABORT(
"MFEM is not built with GSLIB.");
599 size.ProjectCoefficient(size_coeff);
604 size.ProjectCoefficient(size_coeff);
617 disc.ProjectCoefficient(mat_coeff);
624 #ifdef MFEM_USE_GSLIB 627 MFEM_ABORT(
"MFEM is not built with GSLIB.");
635 disc.GetDerivative(1,0,d_x);
636 disc.GetDerivative(1,1,d_y);
639 for (
int i = 0; i < size.Size(); i++)
641 size(i) = std::pow(d_x(i),2)+std::pow(d_y(i),2);
643 const double max = size.Max();
645 for (
int i = 0; i < d_x.Size(); i++)
647 d_x(i) = std::abs(d_x(i));
648 d_y(i) = std::abs(d_y(i));
650 const double eps = 0.01;
651 const double aspr_ratio = 20.0;
652 const double size_ratio = 40.0;
654 for (
int i = 0; i < size.Size(); i++)
656 size(i) = (size(i)/max);
657 aspr(i) = (d_x(i)+eps)/(d_y(i)+eps);
658 aspr(i) = 0.1 + 0.9*(1-size(i))*(1-size(i));
659 if (aspr(i) > aspr_ratio) {aspr(i) = aspr_ratio;}
660 if (aspr(i) < 1.0/aspr_ratio) {aspr(i) = 1.0/aspr_ratio;}
663 const int NE = mesh->
GetNE();
664 double volume = 0.0, volume_ind = 0.0;
666 for (
int i = 0; i < NE; i++)
671 size.GetValues(i, ir, vals);
681 const double avg_zone_size = volume / NE;
683 const double small_avg_ratio = (volume_ind + (volume - volume_ind) /
687 const double small_zone_size = small_avg_ratio * avg_zone_size;
688 const double big_zone_size = size_ratio * small_zone_size;
690 for (
int i = 0; i < size.Size(); i++)
692 const double val = size(i);
693 const double a = (big_zone_size - small_zone_size) / small_zone_size;
694 size(i) = big_zone_size / (1.0+
a*val);
715 #ifdef MFEM_USE_GSLIB 718 MFEM_ABORT(
"MFEM is not built with GSLIB.");
738 #ifdef MFEM_USE_GSLIB 741 MFEM_ABORT(
"MFEM is not built with GSLIB.");
746 size.ProjectCoefficient(size_coeff);
767 default: cout <<
"Unknown target_id: " << target_id << endl;
return 3;
769 if (target_c == NULL)
777 if (metric_combo && bal_expl_combo)
780 metric_combo->ComputeBalancedWeights(x, *target_c, bal_weights);
781 metric_combo->SetWeights(bal_weights);
789 if (barrier_type > 0 || worst_case_type > 0)
797 MFEM_VERIFY(pa ==
false,
"PA for finite differences is not implemented.");
809 default: cout <<
"Unknown quad_type: " << quad_type << endl;
return 3;
814 cout <<
"Triangle quadrature points: " 816 <<
"\nQuadrilateral quadrature points: " 821 cout <<
"Tetrahedron quadrature points: " 823 <<
"\nHexahedron quadrature points: " 825 <<
"\nPrism quadrature points: " 835 if (normalization) { dist = small_phys_size; }
837 if (lim_const != 0.0) { tmop_integ->
EnableLimiting(x0, dist, lim_coeff); }
843 if (adapt_lim_const > 0.0)
845 MFEM_VERIFY(pa ==
false,
"PA is not implemented for adaptive limiting");
850 if (adapt_eval == 0) { adapt_lim_eval =
new AdvectorCG(al); }
851 else if (adapt_eval == 1)
853 #ifdef MFEM_USE_GSLIB 856 MFEM_ABORT(
"MFEM is not built with GSLIB support!");
859 else { MFEM_ABORT(
"Bad interpolation option."); }
882 if (surface_fit_const > 0.0)
884 MFEM_VERIFY(hradaptivity ==
false,
885 "Surface fitting with HR is not implemented yet.");
886 MFEM_VERIFY(pa ==
false,
887 "Surface fitting with PA is not implemented yet.");
892 for (
int i = 0; i < mesh->
GetNE(); i++)
900 for (
int j = 0; j < surf_fit_marker.Size(); j++)
902 if (surf_fit_mat_gf(j) > 0.1 && surf_fit_mat_gf(j) < 0.9)
904 surf_fit_marker[j] =
true;
905 surf_fit_mat_gf(j) = 1.0;
909 surf_fit_marker[j] =
false;
910 surf_fit_mat_gf(j) = 0.0;
914 if (adapt_eval == 0) { adapt_surface =
new AdvectorCG; }
915 else if (adapt_eval == 1)
917 #ifdef MFEM_USE_GSLIB 920 MFEM_ABORT(
"MFEM is not built with GSLIB support!");
923 else { MFEM_ABORT(
"Bad interpolation option."); }
926 surf_fit_coeff, *adapt_surface);
977 else { tmop_integ2 =
new TMOP_Integrator(metric2, target_c, h_metric); }
986 if (lim_const != 0.0) { combo->
EnableLimiting(x0, dist, lim_coeff); }
988 a.AddDomainIntegrator(combo);
992 a.AddDomainIntegrator(tmop_integ);
995 if (pa) {
a.Setup(); }
999 const int NE = mesh->
GetNE();
1000 for (
int i = 0; i < NE; i++)
1008 min_detJ = min(min_detJ, transf->
Jacobian().
Det());
1011 cout <<
"Minimum det(J) of the original mesh is " << min_detJ << endl;
1013 if (min_detJ < 0.0 && barrier_type == 0
1014 && metric_id != 22 && metric_id != 211 && metric_id != 252
1015 && metric_id != 311 && metric_id != 313 && metric_id != 352)
1017 MFEM_ABORT(
"The input mesh is inverted! Try an untangling metric.");
1022 "Untangling is supported only for ideal targets.");
1026 min_detJ /= Wideal.
Det();
1029 min_detJ -= 0.01 * h0.Min();
1033 const double init_energy =
a.GetGridFunctionEnergy(x) /
1034 (hradaptivity ? mesh->
GetNE() : 1);
1035 double init_metric_energy = init_energy;
1036 if (lim_const > 0.0 || adapt_lim_const > 0.0 || surface_fit_const > 0.0)
1040 surf_fit_coeff.constant = 0.0;
1041 init_metric_energy =
a.GetGridFunctionEnergy(x) /
1042 (hradaptivity ? mesh->
GetNE() : 1);
1044 adapt_lim_coeff.
constant = adapt_lim_const;
1045 surf_fit_coeff.constant = surface_fit_const;
1052 char title[] =
"Initial metric values";
1061 if (move_bnd ==
false)
1065 a.SetEssentialBC(ess_bdr);
1070 for (
int i = 0; i < mesh->
GetNBE(); i++)
1074 MFEM_VERIFY(!(
dim == 2 && attr == 3),
1075 "Boundary attribute 3 must be used only for 3D meshes. " 1076 "Adjust the attributes (1/2/3/4 for fixed x/y/z/all " 1077 "components, rest for free nodes), or use -fix-bnd.");
1078 if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
1079 if (attr == 4) { n += nd *
dim; }
1083 for (
int i = 0; i < mesh->
GetNBE(); i++)
1090 for (
int j = 0; j < nd; j++)
1091 { ess_vdofs[n++] = vdofs[j]; }
1095 for (
int j = 0; j < nd; j++)
1096 { ess_vdofs[n++] = vdofs[j+nd]; }
1100 for (
int j = 0; j < nd; j++)
1101 { ess_vdofs[n++] = vdofs[j+2*nd]; }
1105 for (
int j = 0; j < vdofs.
Size(); j++)
1106 { ess_vdofs[n++] = vdofs[j]; }
1109 a.SetEssentialVDofs(ess_vdofs);
1114 Solver *S = NULL, *S_prec = NULL;
1115 const double linsol_rtol = 1e-12;
1116 if (lin_solver == 0)
1118 S =
new DSmoother(1, 1.0, max_lin_iter);
1120 else if (lin_solver == 1)
1137 if (lin_solver == 3 || lin_solver == 4)
1141 MFEM_VERIFY(lin_solver != 4,
"PA l1-Jacobi is not implemented");
1148 auto ds =
new DSmoother((lin_solver == 3) ? 0 : 1, 1.0, 1);
1149 ds->SetPositiveDiagonal(
true);
1162 if (surface_fit_threshold > 0)
1168 if (solver_type == 0)
1178 if (solver_art_type > 0)
1192 x, move_bnd, hradaptivity,
1193 mesh_poly_deg, h_metric_id,
1194 n_hr_iter, n_h_iter);
1196 if (adapt_lim_const > 0.)
1206 ofstream mesh_ofs(
"optimized.mesh");
1207 mesh_ofs.precision(14);
1208 mesh->
Print(mesh_ofs);
1212 const double fin_energy =
a.GetGridFunctionEnergy(x) /
1213 (hradaptivity ? mesh->
GetNE() : 1);
1214 double fin_metric_energy = fin_energy;
1215 if (lim_const > 0.0 || adapt_lim_const > 0.0)
1219 surf_fit_coeff.constant = 0.0;
1220 fin_metric_energy =
a.GetGridFunctionEnergy(x) /
1221 (hradaptivity ? mesh->
GetNE() : 1);
1223 adapt_lim_coeff.
constant = adapt_lim_const;
1224 surf_fit_coeff.constant = surface_fit_const;
1226 std::cout << std::scientific << std::setprecision(4);
1227 cout <<
"Initial strain energy: " << init_energy
1228 <<
" = metrics: " << init_metric_energy
1229 <<
" + extra terms: " << init_energy - init_metric_energy << endl;
1230 cout <<
" Final strain energy: " << fin_energy
1231 <<
" = metrics: " << fin_metric_energy
1232 <<
" + extra terms: " << fin_energy - fin_metric_energy << endl;
1233 cout <<
"The strain energy decreased by: " 1234 << (init_energy - fin_energy) * 100.0 / init_energy <<
" %." << endl;
1239 char title[] =
"Final metric values";
1243 if (adapt_lim_const > 0.0 && visualization)
1247 600, 600, 300, 300);
1251 if (surface_fit_const > 0.0)
1257 600, 900, 300, 300);
1259 900, 900, 300, 300);
1261 double err_avg, err_max;
1263 std::cout <<
"Avg fitting error: " << err_avg << std::endl
1264 <<
"Max fitting error: " << err_max << std::endl;
1271 sock <<
"solution\n";
1276 sock <<
"window_title 'Displacements'\n" 1277 <<
"window_geometry " 1278 << 1200 <<
" " << 0 <<
" " << 600 <<
" " << 600 <<
"\n" 1279 <<
"keys jRmclA" << endl;
1286 delete metric_coeff1;
1287 delete adapt_lim_eval;
1288 delete adapt_surface;
1290 delete hr_adapt_coeff;
1294 delete untangler_metric;
void ComputeUntangleMetricQuantiles(const Vector &x, const FiniteElementSpace &fes)
void discrete_aspr_3d(const Vector &x, Vector &v)
Conjugate gradient method.
int GetNPoints() const
Returns the number of the points in the integration rule.
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
Class for an integration rule - an Array of IntegrationPoint.
3D non-barrier Shape (S) metric.
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
Class for grid function - Vector with associated FE space.
double discrete_size_2d(const Vector &x)
3D barrier Shape+Size (VS) metric, well-posed (invex).
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.
IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto)
void SetPositiveDiagonal(bool pos_diag=true)
Replace diagonal entries with their absolute values.
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
A coefficient that is constant across space and time.
void PrintOptions(std::ostream &out) const
Print the options.
Geometry::Type GetElementBaseGeometry(int i) const
virtual void SetSerialDiscreteTargetOrientation(const GridFunction &tspec_)
int main(int argc, char *argv[])
void EnableNormalization(const GridFunction &x)
Computes the normalization factors of the metric and limiting integrals using the mesh position given...
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
void SetAdaptivityEvaluator(AdaptivityEvaluator *ae)
void PrintUsage(std::ostream &out) const
Print the usage message.
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
int GetAttribute() const
Return element's attribute.
double surface_level_set(const Vector &x)
int Size() const
Returns the size of the vector.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Container class for integration rules.
Data type dense matrix using column-major storage.
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
void GetSurfaceFittingErrors(double &err_avg, double &err_max)
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
int GetNDofs() const
Returns number of degrees of freedom.
void EnableAdaptiveLimiting(const GridFunction &z0, Coefficient &coeff, AdaptivityEvaluator &ae)
Restriction of the node positions to certain regions.
3D non-barrier metric without a type.
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
void EnableNormalization(const GridFunction &x)
Normalization factor that considers all integrators in the combination.
2D barrier Shape+Size (VS) metric (polyconvex).
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
3D barrier Shape+Size (VS) metric, well-posed (invex).
bool Good() const
Return true if the command line options were parsed successfully.
void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, Mesh &mesh, char *title, int position)
void SetExactActionFlag(bool flag_)
Flag to control if exact action of Integration is effected.
void Randomize(int seed=0)
Set random values in the vector.
void SetMinDetPtr(double *md_ptr)
2D barrier Shape+Size (VS) metric (not polyconvex).
3D barrier Shape+Size (VS) metric, well-posed (invex).
int GetNBE() const
Returns number of boundary elements.
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Adds the limiting term to the first integrator. Disables it for the rest.
2D Shifted barrier form of shape metric (mu_2).
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
virtual void SetSerialDiscreteTargetAspectRatio(const GridFunction &tspec_)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
2D barrier shape (S) metric (not polyconvex).
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
void SetTerminationWithMaxSurfaceFittingError(double max_error)
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.
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.
double GetElementSize(ElementTransformation *T, int type=0)
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Limiting of the mesh displacements (general version).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void SetCoefficient(Coefficient &w1)
Sets a scaling Coefficient for the quality metric term of the integrator.
void EnsureNCMesh(bool simplices_nonconforming=false)
3D barrier metric without a type.
Jacobi smoothing for a given bilinear form (no matrix necessary).
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
void SetMaxIter(int max_it)
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Version of QuadraticFECollection with positive basis functions.
double discrete_ori_2d(const Vector &x)
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
void AddFESpaceForUpdate(FiniteElementSpace *fes)
void DiffuseField(ParGridFunction &field, int smooth_steps)
void EnableAdaptiveSurfaceFitting()
2D barrier Shape+Orientation (OS) metric (polyconvex).
A general vector function coefficient.
3D barrier Shape+Size (VS) metric (polyconvex).
void SetAbsTol(double atol)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void AddGridFunctionForUpdate(GridFunction *gf)
void SetRelTol(double rtol)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
double weight_fun(const Vector &x)
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
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)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
void SetAdaptiveLinRtol(const int type=2, const double rtol0=0.5, const double rtol_max=0.9, const double alpha=0.5 *(1.0+sqrt(5.0)), const double gamma=1.0)
Enable adaptive linear solver relative tolerance algorithm.
double adapt_lim_fun(const Vector &x)
int GetDof() const
Returns the number of degrees of freedom in the finite element.
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'.
void AddTMOPIntegrator(TMOP_Integrator *ti)
Adds a new TMOP_Integrator to the combination.
double material_indicator_2d(const Vector &x)
2D barrier Shape+Size (VS) metric (not polyconvex).
int GetNE() const
Returns number of elements.
double discrete_size_3d(const Vector &x)
Class for integration point with weight.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th boundary element.
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
3D Shape (S) metric, untangling version of 303.
IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform)
int GetNBE() const
Returns number of boundary elements in the mesh.
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
int Size() const
Return the logical size of the array.
A general function coefficient.
void SetAttribute(int i, int attr)
Set the attribute of element i.
int DofToVDof(int dof, int vd, int ndofs=-1) const
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
virtual void Print(std::ostream &os=mfem::out) const
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Arbitrary order H1-conforming (continuous) finite elements.
TargetType
Target-matrix construction algorithms supported by this class.
void EnableSurfaceFitting(const GridFunction &s0, const Array< bool > &smarker, Coefficient &coeff, AdaptivityEvaluator &ae)
Fitting of certain DOFs to the zero level set of a function.
const Element * GetBdrElement(int i) const
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
virtual void SetAnalyticTargetSpec(Coefficient *sspec, VectorCoefficient *vspec, TMOPMatrixCoefficient *mspec)
3D barrier Shape+Size (VS) metric (polyconvex).
int material_id(int el_id, const GridFunction &g)
virtual void SetSerialDiscreteTargetSize(const GridFunction &tspec_)
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
2D barrier Shape+Orientation (OS) metric (polyconvex).
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
double GetElementVolume(int i)
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
2D non-barrier metric without a type.
Arbitrary order "L2-conforming" discontinuous finite elements.
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
2D non-barrier Shape+Size+Orientation (VOS) metric (polyconvex).