MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
mesh-optimizer.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12// MFEM Mesh Optimizer Miniapp - Serial/Parallel Shared Code
13
14#include "mfem.hpp"
15#include <fstream>
16#include <iostream>
17
18using namespace mfem;
19using namespace std;
20
22{
23 // semi-circle
24 const real_t xc = x(0) - 0.0, yc = x(1) - 0.5,
25 zc = (x.Size() == 3) ? x(2) - 0.5 : 0.0;
26 const real_t r = sqrt(xc*xc + yc*yc + zc*zc);
27 real_t r1 = 0.45; real_t r2 = 0.55; real_t sf=30.0;
28 real_t val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
29
30 val = std::max((real_t) 0.,val);
31 val = std::min((real_t) 1.,val);
32 return val;
33}
34
35void calc_mass_volume(const GridFunction &g, real_t &mass, real_t &vol)
36{
37 Mesh &mesh = *g.FESpace()->GetMesh();
38 const int NE = mesh.GetNE();
39 Vector g_vals;
40 mass = 0.0, vol = 0.0;
41 for (int e = 0; e < NE; e++)
42 {
45 Tr.OrderJ());
46 g.GetValues(Tr, ir, g_vals);
47 for (int j = 0; j < ir.GetNPoints(); j++)
48 {
49 const IntegrationPoint &ip = ir.IntPoint(j);
50 Tr.SetIntPoint(&ip);
51 mass += g_vals(j) * ip.weight * Tr.Weight();
52 vol += ip.weight * Tr.Weight();
53 }
54 }
55
56#ifdef MFEM_USE_MPI
57 auto gp = dynamic_cast<const ParGridFunction *>(&g);
58 if (gp)
59 {
60 MPI_Comm comm = gp->ParFESpace()->GetComm();
61 MPI_Allreduce(MPI_IN_PLACE, &mass, 1, MPITypeMap<real_t>::mpi_type,
62 MPI_SUM, comm);
63 MPI_Allreduce(MPI_IN_PLACE, &vol, 1, MPITypeMap<real_t>::mpi_type,
64 MPI_SUM, comm);
65 }
66#endif
67}
68
70{
71 // Indicator for small (value -> 1) or big (value -> 0) elements.
72 FunctionCoefficient size_ind_coeff(size_indicator);
73 size.ProjectCoefficient(size_ind_coeff);
74
75 // Determine small/big target sizes based on the total number of
76 // elements and the volume occupied by small elements.
77 real_t volume_ind, volume;
78 calc_mass_volume(size, volume_ind, volume);
79 Mesh &mesh = *size.FESpace()->GetMesh();
80 int NE = mesh.GetNE();
81#ifdef MFEM_USE_MPI
82 auto size_p = dynamic_cast<const ParGridFunction *>(&size);
83 if (size_p) { NE = size_p->ParFESpace()->GetParMesh()->GetGlobalNE(); }
84#endif
85 NCMesh *ncmesh = mesh.ncmesh;
86 // For parallel NC meshes, all tasks have all root elements.
87 NE = (ncmesh) ? ncmesh->GetNumRootElements() : NE;
88 const real_t size_ratio = (mesh.Dimension() == 2) ? 9 : 27;
89 const real_t small_el_size = volume_ind / NE +
90 (volume - volume_ind) / (size_ratio * NE);
91 const real_t big_el_size = size_ratio * small_el_size;
92 for (int i = 0; i < size.Size(); i++)
93 {
94 size(i) = size(i) * small_el_size + (1.0 - size(i)) * big_el_size;
95 }
96}
97
99{
100 real_t xc = x(0)-0.5, yc = x(1)-0.5;
101 real_t th = 22.5*M_PI/180.;
102 real_t xn = cos(th)*xc + sin(th)*yc;
103 real_t yn = -sin(th)*xc + cos(th)*yc;
104 real_t th2 = (th > 45.*M_PI/180) ? M_PI/2 - th : th;
105 real_t stretch = 1/cos(th2);
106 xc = xn/stretch; yc = yn/stretch;
107 real_t tfac = 20;
108 real_t s1 = 3;
109 real_t s2 = 3;
110 real_t wgt = std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) + 1);
111 if (wgt > 1) { wgt = 1; }
112 if (wgt < 0) { wgt = 0; }
113 return wgt;
114}
115
117{
118 return M_PI * x(1) * (1.0 - x(1)) * cos(2 * M_PI * x(0));
119}
120
122{
123 int dim = x.Size();
124 v.SetSize(dim);
125 real_t l1, l2, l3;
126 l1 = 1.;
127 l2 = 1. + 5*x(1);
128 l3 = 1. + 10*x(2);
129 v[0] = l1/pow(l2*l3,0.5);
130 v[1] = l2/pow(l1*l3,0.5);
131 v[2] = l3/pow(l2*l1,0.5);
132}
133
135{
136private:
137 int metric;
138
139public:
140 HessianCoefficient(int dim, int metric_id)
141 : TMOPMatrixCoefficient(dim), metric(metric_id) { }
142
144 const IntegrationPoint &ip)
145 {
146 Vector pos(3);
147 T.Transform(ip, pos);
148 if (metric != 14 && metric != 36 && metric != 85)
149 {
150 const real_t xc = pos(0) - 0.5, yc = pos(1) - 0.5;
151 const real_t r = sqrt(xc*xc + yc*yc);
152 real_t r1 = 0.15; real_t r2 = 0.35; real_t sf=30.0;
153 const real_t eps = 0.5;
154
155 const real_t tan1 = std::tanh(sf*(r-r1)),
156 tan2 = std::tanh(sf*(r-r2));
157
158 K(0, 0) = eps + 1.0 * (tan1 - tan2);
159 K(0, 1) = 0.0;
160 K(1, 0) = 0.0;
161 K(1, 1) = 1.0;
162 }
163 else if (metric == 14 || metric == 36) // Size + Alignment
164 {
165 const real_t xc = pos(0), yc = pos(1);
166 real_t theta = M_PI * yc * (1.0 - yc) * cos(2 * M_PI * xc);
167 real_t alpha_bar = 0.1;
168
169 K(0, 0) = cos(theta);
170 K(1, 0) = sin(theta);
171 K(0, 1) = -sin(theta);
172 K(1, 1) = cos(theta);
173
174 K *= alpha_bar;
175 }
176 else if (metric == 85) // Shape + Alignment
177 {
178 Vector x = pos;
179 real_t xc = x(0)-0.5, yc = x(1)-0.5;
180 real_t th = 22.5*M_PI/180.;
181 real_t xn = cos(th)*xc + sin(th)*yc;
182 real_t yn = -sin(th)*xc + cos(th)*yc;
183 xc = xn; yc=yn;
184
185 real_t tfac = 20;
186 real_t s1 = 3;
187 real_t s2 = 2;
188 real_t wgt = std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) + 1)
189 - std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) - 1);
190 if (wgt > 1) { wgt = 1; }
191 if (wgt < 0) { wgt = 0; }
192
193 xc = pos(0), yc = pos(1);
194 real_t theta = M_PI * (yc) * (1.0 - yc) * cos(2 * M_PI * xc);
195
196 K(0, 0) = cos(theta);
197 K(1, 0) = sin(theta);
198 K(0, 1) = -sin(theta);
199 K(1, 1) = cos(theta);
200
201 real_t asp_ratio_tar = 0.1 + 1*(1-wgt)*(1-wgt);
202
203 K(0, 0) *= 1/pow(asp_ratio_tar,0.5);
204 K(1, 0) *= 1/pow(asp_ratio_tar,0.5);
205 K(0, 1) *= pow(asp_ratio_tar,0.5);
206 K(1, 1) *= pow(asp_ratio_tar,0.5);
207 }
208 }
209
211 const IntegrationPoint &ip, int comp)
212 {
213 Vector pos(3);
214 T.Transform(ip, pos);
215 K = 0.;
216 if (metric != 14 && metric != 85)
217 {
218 const real_t xc = pos(0) - 0.5, yc = pos(1) - 0.5;
219 const real_t r = sqrt(xc*xc + yc*yc);
220 real_t r1 = 0.15; real_t r2 = 0.35; real_t sf=30.0;
221
222 const real_t tan1 = std::tanh(sf*(r-r1)),
223 tan2 = std::tanh(sf*(r-r2));
224 real_t tan1d = 0., tan2d = 0.;
225 if (r > 0.001)
226 {
227 tan1d = (1.-tan1*tan1)*(sf)/r,
228 tan2d = (1.-tan2*tan2)*(sf)/r;
229 }
230
231 K(0, 1) = 0.0;
232 K(1, 0) = 0.0;
233 K(1, 1) = 1.0;
234 if (comp == 0) { K(0, 0) = tan1d*xc - tan2d*xc; }
235 else if (comp == 1) { K(0, 0) = tan1d*yc - tan2d*yc; }
236 }
237 }
238};
239
241{
242private:
243 int dim;
244 // 0 - size target in an annular region,
245 // 1 - size+aspect-ratio in an annular region,
246 // 2 - size+aspect-ratio target for a rotate sine wave.
247 int hr_target_type;
248
249public:
250 HRHessianCoefficient(int dim_, int hr_target_type_ = 0)
251 : TMOPMatrixCoefficient(dim_), dim(dim_),
252 hr_target_type(hr_target_type_) { }
253
255 const IntegrationPoint &ip)
256 {
257 Vector pos(3);
258 T.Transform(ip, pos);
259 if (hr_target_type == 0) // size only circle
260 {
261 real_t small = 0.001, big = 0.01;
262 if (dim == 3) { small = 0.005, big = 0.1; }
263 const real_t xc = pos(0) - 0.5, yc = pos(1) - 0.5;
264 real_t r;
265 if (dim == 2)
266 {
267 r = sqrt(xc*xc + yc*yc);
268 }
269 else
270 {
271 const real_t zc = pos(2) - 0.5;
272 r = sqrt(xc*xc + yc*yc + zc*zc);
273 }
274 real_t r1 = 0.15; real_t r2 = 0.35; real_t sf=30.0;
275
276 const real_t tan1 = std::tanh(sf*(r-r1)),
277 tan2 = std::tanh(sf*(r-r2));
278
279 real_t ind = (tan1 - tan2);
280 if (ind > 1.0) {ind = 1.;}
281 if (ind < 0.0) {ind = 0.;}
282 real_t val = ind * small + (1.0 - ind) * big;
283 K = 0.0;
284 K(0, 0) = 1.0;
285 K(0, 1) = 0.0;
286 K(1, 0) = 0.0;
287 K(1, 1) = 1.0;
288 K(0, 0) *= pow(val,0.5);
289 K(1, 1) *= pow(val,0.5);
290 if (dim == 3) { K(2, 2) = pow(val,0.5); }
291 }
292 else if (hr_target_type == 1) // circle with size and AR
293 {
294 const real_t small = 0.001, big = 0.01;
295 const real_t xc = pos(0)-0.5, yc = pos(1)-0.5;
296 const real_t rv = xc*xc + yc*yc;
297 real_t r = 0;
298 if (rv>0.) {r = sqrt(rv);}
299
300 real_t r1 = 0.2; real_t r2 = 0.3; real_t sf=30.0;
301 const real_t szfac = 1;
302 const real_t asfac = 4;
303 const real_t eps2 = szfac/asfac;
304 const real_t eps1 = szfac;
305
306 real_t tan1 = std::tanh(sf*(r-r1)+1),
307 tan2 = std::tanh(sf*(r-r2)-1);
308 real_t wgt = 0.5*(tan1-tan2);
309
310 tan1 = std::tanh(sf*(r-r1)),
311 tan2 = std::tanh(sf*(r-r2));
312
313 real_t ind = (tan1 - tan2);
314 if (ind > 1.0) {ind = 1.;}
315 if (ind < 0.0) {ind = 0.;}
316 real_t szval = ind * small + (1.0 - ind) * big;
317
318 real_t th = std::atan2(yc,xc)*180./M_PI;
319 if (wgt > 1) { wgt = 1; }
320 if (wgt < 0) { wgt = 0; }
321
322 real_t maxval = eps2 + eps1*(1-wgt)*(1-wgt);
323 real_t minval = eps1;
324 real_t avgval = 0.5*(maxval+minval);
325 real_t ampval = 0.5*(maxval-minval);
326 real_t val1 = avgval + ampval*sin(2.*th*M_PI/180.+90*M_PI/180.);
327 real_t val2 = avgval + ampval*sin(2.*th*M_PI/180.-90*M_PI/180.);
328
329 K(0,1) = 0.0;
330 K(1,0) = 0.0;
331 K(0,0) = val1;
332 K(1,1) = val2;
333
334 K(0,0) *= pow(szval,0.5);
335 K(1,1) *= pow(szval,0.5);
336 }
337 else if (hr_target_type == 2) // sharp rotated sine wave
338 {
339 real_t xc = pos(0)-0.5, yc = pos(1)-0.5;
340 real_t th = 15.5*M_PI/180.;
341 real_t xn = cos(th)*xc + sin(th)*yc;
342 real_t yn = -sin(th)*xc + cos(th)*yc;
343 real_t th2 = (th > 45.*M_PI/180) ? M_PI/2 - th : th;
344 real_t stretch = 1/cos(th2);
345 xc = xn/stretch;
346 yc = yn;
347 real_t tfac = 20;
348 real_t s1 = 3;
349 real_t s2 = 2;
350 real_t yl1 = -0.025;
351 real_t yl2 = 0.025;
352 real_t wgt = std::tanh((tfac*(yc-yl1) + s2*std::sin(s1*M_PI*xc)) + 1) -
353 std::tanh((tfac*(yc-yl2) + s2*std::sin(s1*M_PI*xc)) - 1);
354 if (wgt > 1) { wgt = 1; }
355 if (wgt < 0) { wgt = 0; }
356
357 const real_t eps2 = 25;
358 const real_t eps1 = 1;
359 K(1,1) = eps1/eps2 + eps1*(1-wgt)*(1-wgt);
360 K(0,0) = eps1;
361 K(0,1) = 0.0;
362 K(1,0) = 0.0;
363 }
364 else { MFEM_ABORT("Unsupported option / wrong input."); }
365 }
366
368 const IntegrationPoint &ip, int comp)
369 {
370 K = 0.;
371 }
372};
373
374// Additional IntegrationRules that can be used with the --quad-type option.
377
378// Defined with respect to the icf mesh.
380{
381 const real_t r = sqrt(x(0)*x(0) + x(1)*x(1) + 1e-12);
382 const real_t den = 0.002;
383 real_t l2 = 0.2 + 0.5*std::tanh((r-0.16)/den) - 0.5*std::tanh((r-0.17)/den)
384 + 0.5*std::tanh((r-0.23)/den) - 0.5*std::tanh((r-0.24)/den);
385 return l2;
386}
387
388// Used for the adaptive limiting examples.
390{
391 const real_t xc = x(0) - 0.1, yc = x(1) - 0.2;
392 const real_t r = sqrt(xc*xc + yc*yc);
393 real_t r1 = 0.45; real_t r2 = 0.55; real_t sf=30.0;
394 real_t val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
395
396 val = std::max((real_t) 0.,val);
397 val = std::min((real_t) 1.,val);
398 return val;
399}
400
401// Used for exact surface alignment
403{
404 const int type = 1;
405
406 const int dim = x.Size();
407 if (type == 0)
408 {
409 const real_t sine = 0.25 * std::sin(4 * M_PI * x(0));
410 return (x(1) >= sine + 0.5) ? 1.0 : -1.0;
411 }
412 else
413 {
414 if (dim == 2)
415 {
416 const real_t xc = x(0) - 0.5, yc = x(1) - 0.5;
417 const real_t r = sqrt(xc*xc + yc*yc);
418 return r-0.3;
419 }
420 else
421 {
422 const real_t xc = x(0) - 0.5, yc = x(1) - 0.5, zc = x(2) - 0.5;
423 const real_t r = sqrt(xc*xc + yc*yc + zc*zc);
424 return r-0.3;
425 }
426 }
427}
428
429int material_id(int el_id, const GridFunction &g)
430{
431 const FiniteElementSpace *fes = g.FESpace();
432 const FiniteElement *fe = fes->GetFE(el_id);
433 Vector g_vals;
434 const IntegrationRule &ir =
435 IntRules.Get(fe->GetGeomType(), fes->GetOrder(el_id) + 2);
436
437 real_t integral = 0.0;
438 g.GetValues(el_id, ir, g_vals);
440 int approach = 1;
441 if (approach == 0) // integral based
442 {
443 for (int q = 0; q < ir.GetNPoints(); q++)
444 {
445 const IntegrationPoint &ip = ir.IntPoint(q);
446 Tr->SetIntPoint(&ip);
447 integral += ip.weight * g_vals(q) * Tr->Weight();
448 }
449 return (integral > 0.0) ? 1.0 : 0.0;
450 }
451 else if (approach == 1) // minimum value based
452 {
453 real_t minval = g_vals.Min();
454 return minval > 0.0 ? 1.0 : 0.0;
455 }
456 return 0.0;
457}
458
459void DiffuseField(GridFunction &field, int smooth_steps)
460{
461 // Setup the Laplacian operator
462 BilinearForm *Lap = new BilinearForm(field.FESpace());
464 Lap->Assemble();
465 Lap->Finalize();
466
467 // Setup the smoothing operator
468 DSmoother *S = new DSmoother(0,1.0,smooth_steps);
469 S->iterative_mode = true;
470 S->SetOperator(Lap->SpMat());
471
472 Vector tmp(field.Size());
473 tmp = 0.0;
474 S->Mult(tmp, field);
475
476 delete S;
477 delete Lap;
478}
479
480#ifdef MFEM_USE_MPI
481void DiffuseField(ParGridFunction &field, int smooth_steps)
482{
483 // Setup the Laplacian operator
484 ParBilinearForm *Lap = new ParBilinearForm(field.ParFESpace());
486 Lap->Assemble();
487 Lap->Finalize();
489
490 HypreSmoother *S = new HypreSmoother(*A,0,smooth_steps);
491 S->iterative_mode = true;
492
493 Vector tmp(A->Width());
494 field.SetTrueVector();
495 Vector fieldtrue = field.GetTrueVector();
496 tmp = 0.0;
497 S->Mult(tmp, fieldtrue);
498
499 field.SetFromTrueDofs(fieldtrue);
500
501 delete S;
502 delete A;
503 delete Lap;
504}
505#endif
virtual void EvalGrad(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip, int comp)
Evaluate the derivative of the matrix coefficient with respect to comp in the element described by T ...
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
HRHessianCoefficient(int dim_, int hr_target_type_=0)
virtual void EvalGrad(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip, int comp)
Evaluate the derivative of the matrix coefficient with respect to comp in the element described by T ...
HessianCoefficient(int dim, int metric_id)
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
Data type for scaled Jacobi-type smoother of sparse matrix.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication with Jacobi smoother.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition: eltrans.hpp:131
virtual int OrderJ() const =0
Return the order of the elements of the Jacobian of the transformation.
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:220
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition: fespace.cpp:3168
int GetOrder(int i) const
Returns the polynomial degree of the i'th finite element.
Definition: fespace.hpp:696
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:559
Abstract class for all finite elements.
Definition: fe_base.hpp:239
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:326
A general function coefficient.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:521
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
Definition: gridfunc.hpp:130
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:696
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2346
Wrapper for hypre's ParCSR matrix class.
Definition: hypre.hpp:388
Parallel smoothers in hypre.
Definition: hypre.hpp:1020
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
Definition: hypre.cpp:3777
Class for integration point with weight.
Definition: intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:259
Container class for integration rules.
Definition: intrules.hpp:416
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:1006
Mesh data type.
Definition: mesh.hpp:56
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
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...
Definition: mesh.cpp:357
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:1255
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition: mesh.hpp:291
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1385
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
Definition: ncmesh.hpp:123
int GetNumRootElements()
Return the number of root elements.
Definition: ncmesh.hpp:429
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
Class for parallel bilinear form.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
Class for parallel grid function.
Definition: pgridfunc.hpp:33
void SetFromTrueDofs(const Vector &tv) override
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
@ ClosedUniform
aka closed Newton-Cotes
Definition: intrules.hpp:402
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:686
virtual void SetOperator(const Operator &a)
Set/update the solver for the given operator.
Vector data type.
Definition: vector.hpp:80
int Size() const
Returns the size of the vector.
Definition: vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:538
real_t Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1212
int dim
Definition: ex24.cpp:53
real_t size_indicator(const Vector &x)
int material_id(int el_id, const GridFunction &g)
real_t surface_level_set(const Vector &x)
real_t adapt_lim_fun(const Vector &x)
void DiffuseField(GridFunction &field, int smooth_steps)
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)
void calc_mass_volume(const GridFunction &g, real_t &mass, real_t &vol)
real_t discrete_ori_2d(const Vector &x)
void discrete_aspr_3d(const Vector &x, Vector &v)
float real_t
Definition: config.hpp:43
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:486
Helper struct to convert a C++ type to an MPI type.