119int main (
int argc,
char *argv[])
127 const char *mesh_file =
"icf.mesh";
128 int mesh_poly_deg = 1;
135 real_t adapt_lim_const = 0.0;
139 int solver_iter = 20;
140#ifdef MFEM_USE_SINGLE
141 real_t solver_rtol = 1e-4;
143 real_t solver_rtol = 1e-10;
145 int solver_art_type = 0;
147 int max_lin_iter = 100;
148 bool move_bnd =
true;
150 bool bal_expl_combo =
false;
151 bool hradaptivity =
false;
152 int h_metric_id = -1;
153 bool normalization =
false;
154 bool visualization =
true;
155 int verbosity_level = 0;
156 bool fdscheme =
false;
158 bool exactaction =
false;
159 bool integ_over_targ =
true;
160 const char *devopt =
"cpu";
164 int mesh_node_ordering = 0;
165 int barrier_type = 0;
166 int worst_case_type = 0;
170 args.
AddOption(&mesh_file,
"-m",
"--mesh",
171 "Mesh file to use.");
172 args.
AddOption(&mesh_poly_deg,
"-o",
"--order",
173 "Polynomial degree of mesh finite element space.");
174 args.
AddOption(&rs_levels,
"-rs",
"--refine-serial",
175 "Number of times to refine the mesh uniformly in serial.");
176 args.
AddOption(&rp_levels,
"-rp",
"--refine-parallel",
177 "Number of times to refine the mesh uniformly in parallel.");
178 args.
AddOption(&jitter,
"-ji",
"--jitter",
179 "Random perturbation scaling factor.");
180 args.
AddOption(&metric_id,
"-mid",
"--metric-id",
181 "Mesh optimization metric:\n\t"
183 "1 : |T|^2 -- 2D no type\n\t"
184 "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
185 "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
186 "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
187 "14 : |T-I|^2 -- 2D shape+size+orientation\n\t"
188 "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
189 "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
190 "55 : (tau-1)^2 -- 2D size\n\t"
191 "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
192 "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
193 "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
194 "80 : (1-gamma)mu_2 + gamma mu_77 -- 2D shape+size\n\t"
195 "85 : |T-|T|/sqrt(2)I|^2 -- 2D shape+orientation\n\t"
196 "90 : balanced combo mu_50 & mu_77 -- 2D shape+size\n\t"
197 "94 : balanced combo mu_2 & mu_56 -- 2D shape+size\n\t"
198 "98 : (1/tau)|T-I|^2 -- 2D shape+size+orientation\n\t"
201 "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
202 "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
203 "303: (|T|^2)/3/tau^(2/3)-1 -- 3D shape\n\t"
204 "304: (|T|^3)/3^{3/2}/tau-1 -- 3D shape\n\t"
206 "313: (|T|^2)(tau-tau0)^(-2/3)/3 -- 3D untangling\n\t"
207 "315: (tau-1)^2 -- 3D no type\n\t"
208 "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D no type\n\t"
209 "321: |T-T^-t|^2 -- 3D shape+size\n\t"
210 "322: |T-adjT^-t|^2 -- 3D shape+size\n\t"
211 "323: |J|^3-3sqrt(3)ln(det(J))-3sqrt(3) -- 3D shape+size\n\t"
212 "328: balanced combo mu_301 & mu_316 -- 3D shape+size\n\t"
213 "332: (1-gamma) mu_302 + gamma mu_315 -- 3D shape+size\n\t"
214 "333: (1-gamma) mu_302 + gamma mu_316 -- 3D shape+size\n\t"
215 "334: (1-gamma) mu_303 + gamma mu_316 -- 3D shape+size\n\t"
216 "328: balanced combo mu_302 & mu_318 -- 3D shape+size\n\t"
217 "347: (1-gamma) mu_304 + gamma mu_316 -- 3D shape+size\n\t"
219 "360: (|T|^3)/3^{3/2}-tau -- 3D shape\n\t"
221 "11 : (1/4*alpha)|A-(adjA)^T(W^TW)/omega|^2 -- 2D shape\n\t"
222 "36 : (1/alpha)|A-W|^2 -- 2D shape+size+orientation\n\t"
223 "107: (1/2*alpha)|A-|A|/|W|W|^2 -- 2D shape+orientation\n\t"
224 "126: (1-gamma)nu_11 + gamma*nu_14a -- 2D shape+size\n\t"
226 args.
AddOption(&target_id,
"-tid",
"--target-id",
227 "Target (ideal element) type:\n\t"
228 "1: Ideal shape, unit size\n\t"
229 "2: Ideal shape, equal size\n\t"
230 "3: Ideal shape, initial size\n\t"
231 "4: Given full analytic Jacobian (in physical space)\n\t"
232 "5: Ideal shape, given size (in physical space)");
233 args.
AddOption(&lim_const,
"-lc",
"--limit-const",
"Limiting constant.");
234 args.
AddOption(&adapt_lim_const,
"-alc",
"--adapt-limit-const",
235 "Adaptive limiting coefficient constant.");
236 args.
AddOption(&quad_type,
"-qt",
"--quad-type",
237 "Quadrature rule type:\n\t"
238 "1: Gauss-Lobatto\n\t"
239 "2: Gauss-Legendre\n\t"
240 "3: Closed uniform points");
241 args.
AddOption(&quad_order,
"-qo",
"--quad_order",
242 "Order of the quadrature rule.");
243 args.
AddOption(&solver_type,
"-st",
"--solver-type",
244 " Type of solver: (default) 0: Newton, 1: LBFGS");
245 args.
AddOption(&solver_iter,
"-ni",
"--newton-iters",
246 "Maximum number of Newton iterations.");
247 args.
AddOption(&solver_rtol,
"-rtol",
"--newton-rel-tolerance",
248 "Relative tolerance for the Newton solver.");
249 args.
AddOption(&solver_art_type,
"-art",
"--adaptive-rel-tol",
250 "Type of adaptive relative linear solver tolerance:\n\t"
251 "0: None (default)\n\t"
252 "1: Eisenstat-Walker type 1\n\t"
253 "2: Eisenstat-Walker type 2");
254 args.
AddOption(&lin_solver,
"-ls",
"--lin-solver",
259 "3: MINRES + Jacobi preconditioner\n\t"
260 "4: MINRES + l1-Jacobi preconditioner");
261 args.
AddOption(&max_lin_iter,
"-li",
"--lin-iter",
262 "Maximum number of iterations in the linear solve.");
263 args.
AddOption(&move_bnd,
"-bnd",
"--move-boundary",
"-fix-bnd",
265 "Enable motion along horizontal and vertical boundaries.");
266 args.
AddOption(&combomet,
"-cmb",
"--combo-type",
267 "Combination of metrics options:\n\t"
268 "0: Use single metric\n\t"
269 "1: Shape + space-dependent size given analytically\n\t"
270 "2: Shape + adapted size given discretely; shared target");
271 args.
AddOption(&bal_expl_combo,
"-bec",
"--balance-explicit-combo",
272 "-no-bec",
"--balance-explicit-combo",
273 "Automatic balancing of explicit combo metrics.");
274 args.
AddOption(&hradaptivity,
"-hr",
"--hr-adaptivity",
"-no-hr",
275 "--no-hr-adaptivity",
276 "Enable hr-adaptivity.");
277 args.
AddOption(&h_metric_id,
"-hmid",
"--h-metric",
278 "Same options as metric_id. Used to determine refinement"
279 " type for each element if h-adaptivity is enabled.");
280 args.
AddOption(&normalization,
"-nor",
"--normalization",
"-no-nor",
281 "--no-normalization",
282 "Make all terms in the optimization functional unitless.");
283 args.
AddOption(&fdscheme,
"-fd",
"--fd_approximation",
284 "-no-fd",
"--no-fd-approx",
285 "Enable finite difference based derivative computations.");
286 args.
AddOption(&exactaction,
"-ex",
"--exact_action",
287 "-no-ex",
"--no-exact-action",
288 "Enable exact action of TMOP_Integrator.");
289 args.
AddOption(&integ_over_targ,
"-it",
"--integrate-target",
290 "-ir",
"--integrate-reference",
291 "Integrate over target (-it) or reference (-ir) element.");
292 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
293 "--no-visualization",
294 "Enable or disable GLVis visualization.");
295 args.
AddOption(&verbosity_level,
"-vl",
"--verbosity-level",
296 "Verbosity level for the involved iterative solvers:\n\t"
298 "1: Newton iterations\n\t"
299 "2: Newton iterations + linear solver summaries\n\t"
300 "3: newton iterations + linear solver iterations");
301 args.
AddOption(&adapt_eval,
"-ae",
"--adaptivity-evaluator",
302 "0 - Advection based (DEFAULT), 1 - GSLIB.");
303 args.
AddOption(&devopt,
"-d",
"--device",
304 "Device configuration string, see Device::Configure().");
305 args.
AddOption(&pa,
"-pa",
"--partial-assembly",
"-no-pa",
306 "--no-partial-assembly",
"Enable Partial Assembly.");
307 args.
AddOption(&n_hr_iter,
"-nhr",
"--n_hr_iter",
308 "Number of hr-adaptivity iterations.");
309 args.
AddOption(&n_h_iter,
"-nh",
"--n_h_iter",
310 "Number of h-adaptivity iterations per r-adaptivity"
312 args.
AddOption(&mesh_node_ordering,
"-mno",
"--mesh_node_ordering",
313 "Ordering of mesh nodes."
314 "0 (default): byNodes, 1: byVDIM");
315 args.
AddOption(&barrier_type,
"-btype",
"--barrier-type",
317 "1 - Shifted Barrier,"
318 "2 - Pseudo Barrier.");
319 args.
AddOption(&worst_case_type,
"-wctype",
"--worst-case-type",
331 if (h_metric_id < 0) { h_metric_id = metric_id; }
335 MFEM_VERIFY(strcmp(devopt,
"cpu")==0,
"HR-adaptivity is currently only"
336 " supported on cpus.");
339 if (myid == 0) { device.
Print();}
342 Mesh *mesh =
new Mesh(mesh_file, 1, 1,
false);
343 for (
int lev = 0; lev < rs_levels; lev++)
352 for (
int lev = 0; lev < rp_levels; lev++)
362 if (mesh_poly_deg <= 0)
391 for (
int i = 0; i < pmesh->
GetNE(); i++)
397 for (
int j = 0; j < dofs.
Size(); j++)
399 h0(dofs[j]) = min(h0(dofs[j]), hi);
405 MPI_SUM, MPI_COMM_WORLD);
406 const real_t small_phys_size = pow(vol_glb, 1.0 /
dim) / 100.0;
419 for (
int i = 0; i < pfespace->
GetNDofs(); i++)
421 for (
int d = 0; d <
dim; d++)
427 for (
int i = 0; i < pfespace->
GetNBE(); i++)
432 for (
int j = 0; j < vdofs.
Size(); j++) { rdm(vdofs[j]) = 0.0; }
443 ostringstream mesh_name;
444 mesh_name <<
"perturbed.mesh";
445 ofstream mesh_ofs(mesh_name.str().c_str());
446 mesh_ofs.precision(8);
505 if (myid == 0) { cout <<
"Unknown metric_id: " << metric_id << endl; }
524 default: cout <<
"Metric_id not supported for h-adaptivity: " << h_metric_id <<
531 switch (barrier_type)
533 case 0: btype = TMOP_WorstCaseUntangleOptimizer_Metric::BarrierType::None;
535 case 1: btype = TMOP_WorstCaseUntangleOptimizer_Metric::BarrierType::Shifted;
537 case 2: btype = TMOP_WorstCaseUntangleOptimizer_Metric::BarrierType::Pseudo;
539 default: cout <<
"barrier_type not supported: " << barrier_type << endl;
544 switch (worst_case_type)
546 case 0: wctype = TMOP_WorstCaseUntangleOptimizer_Metric::WorstCaseType::None;
548 case 1: wctype = TMOP_WorstCaseUntangleOptimizer_Metric::WorstCaseType::Beta;
550 case 2: wctype = TMOP_WorstCaseUntangleOptimizer_Metric::WorstCaseType::PMean;
552 default: cout <<
"worst_case_type not supported: " << worst_case_type << endl;
557 if (barrier_type > 0 || worst_case_type > 0)
559 if (barrier_type > 0)
561 MFEM_VERIFY(metric_id == 4 || metric_id == 14 || metric_id == 66,
562 "Metric not supported for shifted/pseudo barriers.");
573 if (metric_id < 300 || h_metric_id < 300)
575 MFEM_VERIFY(
dim == 2,
"Incompatible metric for 3D meshes");
577 if (metric_id >= 300 || h_metric_id >= 300)
579 MFEM_VERIFY(
dim == 3,
"Incompatible metric for 2D meshes");
593 pa ? AssemblyLevel::PARTIAL : AssemblyLevel::LEGACY;
622 MFEM_ABORT(
"MFEM is not built with GSLIB.");
637 disc.ProjectCoefficient(mat_coeff);
647 MFEM_ABORT(
"MFEM is not built with GSLIB.");
654 disc.GetDerivative(1,0,d_x);
655 disc.GetDerivative(1,1,d_y);
658 for (
int i = 0; i < size.Size(); i++)
660 size(i) = std::pow(d_x(i),2)+std::pow(d_y(i),2);
662 const real_t max = size.Max();
665 MPI_MAX, MPI_COMM_WORLD);
667 for (
int i = 0; i < d_x.Size(); i++)
669 d_x(i) = std::abs(d_x(i));
670 d_y(i) = std::abs(d_y(i));
673 const real_t aspr_ratio = 20.0;
674 const real_t size_ratio = 40.0;
676 for (
int i = 0; i < size.Size(); i++)
678 size(i) = (size(i)/max_all);
679 aspr(i) = (d_x(i)+eps)/(d_y(i)+eps);
680 aspr(i) = 0.1 + 0.9*(1-size(i))*(1-size(i));
681 if (aspr(i) > aspr_ratio) {aspr(i) = aspr_ratio;}
682 if (aspr(i) < 1.0/aspr_ratio) {aspr(i) = 1.0/aspr_ratio;}
685 const int NE = pmesh->
GetNE();
686 real_t volume = 0.0, volume_ind = 0.0;
688 for (
int i = 0; i < NE; i++)
693 size.GetValues(i, ir, vals);
702 real_t volume_all, volume_ind_all;
704 MPI_SUM, MPI_COMM_WORLD);
705 MPI_Allreduce(&volume_ind, &volume_ind_all, 1,
709 const real_t avg_zone_size = volume_all / NE_ALL;
711 const real_t small_avg_ratio =
712 (volume_ind_all + (volume_all - volume_ind_all) / size_ratio)
715 const real_t small_zone_size = small_avg_ratio * avg_zone_size;
716 const real_t big_zone_size = size_ratio * small_zone_size;
718 for (
int i = 0; i < size.Size(); i++)
720 const real_t val = size(i);
721 const real_t a = (big_zone_size - small_zone_size) / small_zone_size;
722 size(i) = big_zone_size / (1.0+
a*val);
746 MFEM_ABORT(
"MFEM is not built with GSLIB.");
768 MFEM_ABORT(
"MFEM is not built with GSLIB.");
773 size.ProjectCoefficient(size_coeff);
795 if (myid == 0) { cout <<
"Unknown target_id: " << target_id << endl; }
798 if (target_c == NULL)
806 if (metric_combo && bal_expl_combo)
809 metric_combo->ComputeBalancedWeights(x, *target_c, bal_weights);
810 metric_combo->SetWeights(bal_weights);
816 auto tmop_integ =
new TMOP_Integrator(metric_to_use, target_c, h_metric);
817 tmop_integ->IntegrateOverTarget(integ_over_targ);
818 if (barrier_type > 0 || worst_case_type > 0)
820 tmop_integ->ComputeUntangleMetricQuantiles(x, *pfespace);
826 MFEM_VERIFY(pa ==
false,
"PA for finite differences is not implemented.");
827 tmop_integ->EnableFiniteDifferences(x);
829 tmop_integ->SetExactActionFlag(exactaction);
839 if (myid == 0) { cout <<
"Unknown quad_type: " << quad_type << endl; }
842 tmop_integ->SetIntegrationRules(*irules, quad_order);
843 if (myid == 0 &&
dim == 2)
845 cout <<
"Triangle quadrature points: "
847 <<
"\nQuadrilateral quadrature points: "
850 if (myid == 0 &&
dim == 3)
852 cout <<
"Tetrahedron quadrature points: "
854 <<
"\nHexahedron quadrature points: "
856 <<
"\nPrism quadrature points: "
866 if (normalization) { dist = small_phys_size; }
868 if (lim_const != 0.0) { tmop_integ->EnableLimiting(x0, dist, lim_coeff); }
874 if (adapt_lim_const > 0.0)
876 MFEM_VERIFY(pa ==
false,
"PA is not implemented for adaptive limiting");
881 if (adapt_eval == 0) { adapt_lim_eval =
new AdvectorCG(al); }
882 else if (adapt_eval == 1)
887 MFEM_ABORT(
"MFEM is not built with GSLIB support!");
890 else { MFEM_ABORT(
"Bad interpolation option."); }
892 tmop_integ->EnableAdaptiveLimiting(adapt_lim_gf0, adapt_lim_coeff,
904 if (normalization) { tmop_integ->ParEnableNormalization(x0); }
913 if (pa) {
a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
924 tmop_integ->SetCoefficient(*metric_coeff1);
939 else { tmop_integ2 =
new TMOP_Integrator(metric2, target_c, h_metric); }
949 if (lim_const != 0.0) { combo->
EnableLimiting(x0, dist, lim_coeff); }
951 a.AddDomainIntegrator(combo);
955 a.AddDomainIntegrator(tmop_integ);
958 if (pa) {
a.Setup(); }
962 const int NE = pmesh->
GetNE();
963 for (
int i = 0; i < NE; i++)
971 min_detJ = min(min_detJ, transf->
Jacobian().
Det());
976 MPI_MIN, MPI_COMM_WORLD);
979 { cout <<
"Minimum det(J) of the original mesh is " << min_detJ << endl; }
981 if (min_detJ < 0.0 && barrier_type == 0
982 && metric_id != 22 && metric_id != 211 && metric_id != 252
983 && metric_id != 311 && metric_id != 313 && metric_id != 352)
985 MFEM_ABORT(
"The input mesh is inverted! Try an untangling metric.");
990 "Untangling is supported only for ideal targets.");
994 min_detJ /= Wideal.
Det();
998 MPI_MIN, MPI_COMM_WORLD);
1000 min_detJ -= 0.01 * h0min_all;
1004 const real_t init_energy =
a.GetParGridFunctionEnergy(x) /
1006 real_t init_metric_energy = init_energy;
1007 if (lim_const > 0.0 || adapt_lim_const > 0.0)
1011 init_metric_energy =
a.GetParGridFunctionEnergy(x) /
1014 adapt_lim_coeff.
constant = adapt_lim_const;
1021 char title[] =
"Initial metric values";
1029 if (move_bnd ==
false)
1033 a.SetEssentialBC(ess_bdr);
1038 for (
int i = 0; i < pmesh->
GetNBE(); i++)
1042 MFEM_VERIFY(!(
dim == 2 && attr == 3),
1043 "Boundary attribute 3 must be used only for 3D meshes. "
1044 "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
1045 "components, rest for free nodes), or use -fix-bnd.");
1046 if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
1047 if (attr == 4) { n += nd *
dim; }
1051 for (
int i = 0; i < pmesh->
GetNBE(); i++)
1058 for (
int j = 0; j < nd; j++)
1059 { ess_vdofs[n++] = vdofs[j]; }
1063 for (
int j = 0; j < nd; j++)
1064 { ess_vdofs[n++] = vdofs[j+nd]; }
1068 for (
int j = 0; j < nd; j++)
1069 { ess_vdofs[n++] = vdofs[j+2*nd]; }
1073 for (
int j = 0; j < vdofs.
Size(); j++)
1074 { ess_vdofs[n++] = vdofs[j]; }
1077 a.SetEssentialVDofs(ess_vdofs);
1082 Solver *S = NULL, *S_prec = NULL;
1083#ifdef MFEM_USE_SINGLE
1084 const real_t linsol_rtol = 1e-5;
1086 const real_t linsol_rtol = 1e-12;
1090 if (verbosity_level == 2)
1092 if (verbosity_level > 2)
1094 if (lin_solver == 0)
1096 S =
new DSmoother(1, 1.0, max_lin_iter);
1098 else if (lin_solver == 1)
1114 if (lin_solver == 3 || lin_solver == 4)
1118 MFEM_VERIFY(lin_solver != 4,
"PA l1-Jacobi is not implemented");
1128 hs->SetPositiveDiagonal(
true);
1151 if (solver_art_type > 0)
1157 if (verbosity_level > 0)
1168 x, move_bnd, hradaptivity,
1169 mesh_poly_deg, h_metric_id,
1170 n_hr_iter, n_h_iter);
1172 if (adapt_lim_const > 0.)
1182 ostringstream mesh_name;
1183 mesh_name <<
"optimized.mesh";
1184 ofstream mesh_ofs(mesh_name.str().c_str());
1185 mesh_ofs.precision(8);
1190 const real_t fin_energy =
a.GetParGridFunctionEnergy(x) /
1192 real_t fin_metric_energy = fin_energy;
1193 if (lim_const > 0.0 || adapt_lim_const > 0.0)
1197 fin_metric_energy =
a.GetParGridFunctionEnergy(x) /
1200 adapt_lim_coeff.
constant = adapt_lim_const;
1204 std::cout << std::scientific << std::setprecision(4);
1205 cout <<
"Initial strain energy: " << init_energy
1206 <<
" = metrics: " << init_metric_energy
1207 <<
" + extra terms: " << init_energy - init_metric_energy << endl;
1208 cout <<
" Final strain energy: " << fin_energy
1209 <<
" = metrics: " << fin_metric_energy
1210 <<
" + extra terms: " << fin_energy - fin_metric_energy << endl;
1211 cout <<
"The strain energy decreased by: "
1212 << (init_energy - fin_energy) * 100.0 / init_energy <<
" %." << endl;
1218 char title[] =
"Final metric values";
1222 if (adapt_lim_const > 0.0 && visualization)
1226 600, 600, 300, 300);
1236 sock.
open(
"localhost", 19916);
1237 sock <<
"solution\n";
1243 sock <<
"window_title 'Displacements'\n"
1244 <<
"window_geometry "
1245 << 1200 <<
" " << 0 <<
" " << 600 <<
" " << 600 <<
"\n"
1246 <<
"keys jRmclA" << endl;
1254 delete metric_coeff1;
1255 delete adapt_lim_eval;
1257 delete hr_adapt_coeff;
1261 delete untangler_metric;
virtual void SetAnalyticTargetSpec(Coefficient *sspec, VectorCoefficient *vspec, TMOPMatrixCoefficient *mspec)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
int Size() const
Return the logical size of the array.
Conjugate gradient method.
A coefficient that is constant across space and time.
Data type for scaled Jacobi-type smoother of sparse matrix.
Data type dense matrix using column-major storage.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
virtual void SetParDiscreteTargetAspectRatio(const ParGridFunction &tspec_)
virtual void SetParDiscreteTargetSize(const ParGridFunction &tspec_)
void SetAdaptivityEvaluator(AdaptivityEvaluator *ae)
virtual void SetParDiscreteTargetOrientation(const ParGridFunction &tspec_)
int GetAttribute() const
Return element's attribute.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
int GetNBE() const
Returns number of boundary elements in the mesh.
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
int GetDof() const
Returns the number of degrees of freedom in the finite element.
A general function coefficient.
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Arbitrary order H1-conforming (continuous) finite elements.
Parallel smoothers in hypre.
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
@ l1Jacobi
l1-scaled Jacobi
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
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.
Container class for integration rules.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void SetRelTol(real_t rtol)
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)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
int GetNE() const
Returns number of elements.
int Dimension() const
Dimension of the reference space used within the elements.
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
real_t GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
real_t GetElementVolume(int i)
int GetNBE() const
Returns number of boundary elements.
long long GetGlobalNE() const
Return the total (global) number of elements.
void EnsureNCMesh(bool simplices_nonconforming=false)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Geometry::Type GetElementBaseGeometry(int i) const
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void SetAdaptiveLinRtol(const int type=2, const real_t rtol0=0.5, const real_t rtol_max=0.9, const real_t alpha=0.5 *(1.0+sqrt(5.0)), const real_t gamma=1.0)
Enable adaptive linear solver relative tolerance algorithm.
Jacobi smoothing for a given bilinear form (no matrix necessary).
void SetPositiveDiagonal(bool pos_diag=true)
Replace diagonal entries with their absolute values.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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...
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
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 FiniteElement * GetFE(int i) const override
Class for parallel grid function.
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void SaveAsOne(const char *fname, int precision=16) const
Class for parallel meshes.
void SetNodalFESpace(FiniteElementSpace *nfes) override
void PrintAsOne(std::ostream &out=mfem::out, const std::string &comments="") const
Version of QuadraticFECollection with positive basis functions.
void ParEnableNormalization(const ParGridFunction &x)
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.
void AddTMOPIntegrator(TMOP_Integrator *ti)
Adds a new TMOP_Integrator to the combination.
void AddGridFunctionForUpdate(GridFunction *gf)
void AddFESpaceForUpdate(FiniteElementSpace *fes)
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
void SetMinDetPtr(real_t *md_ptr)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
2D barrier Shape+Orientation (OS) metric (polyconvex).
2D barrier Shape+Size (VS) metric (polyconvex).
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
void SetExactActionFlag(bool flag_)
Flag to control if exact action of Integration is effected.
void SetCoefficient(Coefficient &w1)
Sets a scaling Coefficient for the quality metric term of the integrator.
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
void IntegrateOverTarget(bool integ_over_target_)
2D non-barrier metric without a type.
2D barrier Shape+Size (VS) metric (not polyconvex).
2D barrier Shape+Size (VS) metric (not polyconvex).
2D non-barrier Shape+Size+Orientation (VOS) metric (polyconvex).
2D Shifted barrier form of shape metric (mu_2).
2D barrier shape (S) metric (not polyconvex).
2D barrier Shape+Orientation (OS) metric (polyconvex).
2D compound barrier Shape+Size (VS) metric (balanced).
2D compound barrier Shape+Size (VS) metric (balanced).
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
3D Shape (S) metric, untangling version of 303.
3D barrier Shape+Size (VS) metric, well-posed (invex).
3D barrier Shape+Size (VS) metric, well-posed (invex).
3D barrier Shape+Size (VS) metric, well-posed (invex).
3D compound barrier Shape+Size (VS) metric (polyconvex, balanced).
3D compound barrier Shape+Size (VS) metric (polyconvex).
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
3D compound barrier Shape+Size (VS) metric (polyconvex, balanced).
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
3D non-barrier Shape (S) metric.
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...
void SetVolumeScale(real_t vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
TargetType
Target-matrix construction algorithms supported by this class.
A general vector function coefficient.
void Randomize(int seed=0)
Set random values in the vector.
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
real_t Min() const
Returns the minimal element of the vector.
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
real_t adapt_lim_fun(const Vector &x)
IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform)
real_t weight_fun(const Vector &x)
real_t material_indicator_2d(const Vector &x)
IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto)
void ConstructSizeGF(GridFunction &size)
real_t discrete_ori_2d(const Vector &x)
void discrete_aspr_3d(const Vector &x, Vector &v)
void DiffuseField(ParGridFunction &field, int smooth_steps)
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)
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Settings for the output behavior of the IterativeSolver.
PrintLevel & Iterations()
PrintLevel & FirstAndLast()
Helper struct to convert a C++ type to an MPI type.