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