MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
pfem_extras.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 #include "pfem_extras.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 using namespace std;
17 
18 namespace mfem
19 {
20 
21 namespace common
22 {
23 
24 H1_ParFESpace::H1_ParFESpace(ParMesh *m,
25  const int p, const int space_dim, const int type,
26  int vdim, int order)
27  : ParFiniteElementSpace(m, new H1_FECollection(p,space_dim,type),vdim,order)
28 {
29  FEC_ = this->FiniteElementSpace::fec;
30 }
31 
33 {
34  delete FEC_;
35 }
36 
37 ND_ParFESpace::ND_ParFESpace(ParMesh *m, const int p, const int space_dim,
38  int vdim, int order)
39  : ParFiniteElementSpace(m, new ND_FECollection(p,space_dim),vdim,order)
40 {
41  FEC_ = this->FiniteElementSpace::fec;
42 }
43 
45 {
46  delete FEC_;
47 }
48 
49 RT_ParFESpace::RT_ParFESpace(ParMesh *m, const int p, const int space_dim,
50  int vdim, int order)
51  : ParFiniteElementSpace(m, new RT_FECollection(p-1,space_dim),vdim,order)
52 {
53  FEC_ = this->FiniteElementSpace::fec;
54 }
55 
57 {
58  delete FEC_;
59 }
60 
61 L2_ParFESpace::L2_ParFESpace(ParMesh *m, const int p, const int space_dim,
62  int vdim, int order)
63  : ParFiniteElementSpace(m, new L2_FECollection(p,space_dim),vdim,order)
64 {
65  FEC_ = this->FiniteElementSpace::fec;
66 }
67 
69 {
70  delete FEC_;
71 }
72 
74 {}
75 
79 {
81 }
82 
86 {
88 }
89 
93 {
95 }
96 
99  ParFiniteElementSpace & HCurlFESpace,
100  const int & irOrder,
101  ParBilinearForm * s0,
102  ParMixedBilinearForm * weakDiv,
104  : H1FESpace_(&H1FESpace),
105  HCurlFESpace_(&HCurlFESpace),
106  s0_(s0),
107  weakDiv_(weakDiv),
108  grad_(grad),
109  psi_(NULL),
110  xDiv_(NULL),
111  S0_(NULL),
112  amg_(NULL),
113  pcg_(NULL),
114  ownsS0_(s0 == NULL),
115  ownsWeakDiv_(weakDiv == NULL),
116  ownsGrad_(grad == NULL)
117 {
118  ess_bdr_.SetSize(H1FESpace_->GetParMesh()->bdr_attributes.Max());
119  ess_bdr_ = 1;
120  H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
121 
122  int geom = H1FESpace_->GetFE(0)->GetGeomType();
123  const IntegrationRule * ir = &IntRules.Get(geom, irOrder);
124 
125  if ( s0 == NULL )
126  {
127  s0_ = new ParBilinearForm(H1FESpace_);
129  diffInteg->SetIntRule(ir);
130  s0_->AddDomainIntegrator(diffInteg);
131  s0_->Assemble();
132  s0_->Finalize();
133  S0_ = new HypreParMatrix;
134  }
135  if ( weakDiv_ == NULL )
136  {
137  weakDiv_ = new ParMixedBilinearForm(HCurlFESpace_, H1FESpace_);
139  wdivInteg->SetIntRule(ir);
140  weakDiv_->AddDomainIntegrator(wdivInteg);
141  weakDiv_->Assemble();
142  weakDiv_->Finalize();
143  }
144  if ( grad_ == NULL )
145  {
146  grad_ = new ParDiscreteGradOperator(H1FESpace_, HCurlFESpace_);
147  grad_->Assemble();
148  grad_->Finalize();
149  }
150 
151  psi_ = new ParGridFunction(H1FESpace_);
152  xDiv_ = new ParGridFunction(H1FESpace_);
153 }
154 
156 {
157  delete psi_;
158  delete xDiv_;
159 
160  delete amg_;
161  delete pcg_;
162 
163  delete S0_;
164 
165  delete s0_;
166  delete weakDiv_;
167 }
168 
169 void
170 IrrotationalProjector::InitSolver() const
171 {
172  delete pcg_;
173  delete amg_;
174 
175  amg_ = new HypreBoomerAMG(*S0_);
176  amg_->SetPrintLevel(0);
177  pcg_ = new HyprePCG(*S0_);
178  pcg_->SetTol(1e-14);
179  pcg_->SetMaxIter(200);
180  pcg_->SetPrintLevel(0);
181  pcg_->SetPreconditioner(*amg_);
182 }
183 
184 void
186 {
187  // Compute the divergence of x
188  weakDiv_->Mult(x,*xDiv_); *xDiv_ *= -1.0;
189 
190  // Apply essential BC and form linear system
191  *psi_ = 0.0;
192  s0_->FormLinearSystem(ess_bdr_tdofs_, *psi_, *xDiv_, *S0_, Psi_, RHS_);
193 
194  // Solve the linear system for Psi
195  if ( pcg_ == NULL ) { this->InitSolver(); }
196  pcg_->Mult(RHS_, Psi_);
197 
198  // Compute the parallel grid function corresponding to Psi
199  s0_->RecoverFEMSolution(Psi_, *xDiv_, *psi_);
200 
201  // Compute the irrotational portion of x
202  grad_->Mult(*psi_, y);
203 }
204 
205 void
207 {
208  delete pcg_; pcg_ = NULL;
209  delete amg_; amg_ = NULL;
210  delete S0_; S0_ = new HypreParMatrix;
211 
212  psi_->Update();
213  xDiv_->Update();
214 
215  if ( ownsS0_ )
216  {
217  s0_->Update();
218  s0_->Assemble();
219  s0_->Finalize();
220  }
221  if ( ownsWeakDiv_ )
222  {
223  weakDiv_->Update();
224  weakDiv_->Assemble();
225  weakDiv_->Finalize();
226  }
227  if ( ownsGrad_ )
228  {
229  grad_->Update();
230  grad_->Assemble();
231  grad_->Finalize();
232  }
233 
234  H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
235 }
236 
239  ParFiniteElementSpace & HCurlFESpace,
240  const int & irOrder,
241  ParBilinearForm * s0,
242  ParMixedBilinearForm * weakDiv,
244  : IrrotationalProjector(H1FESpace,HCurlFESpace, irOrder, s0, weakDiv, grad)
245 {}
246 
248 {}
249 
250 void
252 {
253  this->IrrotationalProjector::Mult(x, y);
254  y -= x;
255  y *= -1.0;
256 }
257 
258 void
260 {
262 }
263 
264 void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
265  ParMesh &pmesh, const char *title,
266  int x, int y, int w, int h, const char *keys)
267 {
268  MPI_Comm comm = pmesh.GetComm();
269 
270  int num_procs, myid;
271  MPI_Comm_size(comm, &num_procs);
272  MPI_Comm_rank(comm, &myid);
273 
274  bool newly_opened = false;
275  int connection_failed;
276 
277  do
278  {
279  if (myid == 0)
280  {
281  if (!sock.is_open() || !sock)
282  {
283  sock.open(vishost, visport);
284  sock.precision(8);
285  newly_opened = true;
286  }
287  sock << "mesh\n";
288  }
289 
290  pmesh.PrintAsOne(sock);
291 
292  if (myid == 0 && newly_opened)
293  {
294  sock << "window_title '" << title << "'\n"
295  << "window_geometry "
296  << x << " " << y << " " << w << " " << h << "\n";
297  if ( keys ) { sock << "keys " << keys << "\n"; }
298  sock << endl;
299  }
300 
301  if (myid == 0)
302  {
303  connection_failed = !sock && !newly_opened;
304  }
305  MPI_Bcast(&connection_failed, 1, MPI_INT, 0, comm);
306  }
307  while (connection_failed);
308 }
309 
310 void VisualizeField(socketstream &sock, const char *vishost, int visport,
311  const ParGridFunction &gf, const char *title,
312  int x, int y, int w, int h, const char *keys, bool vec)
313 {
314  ParMesh &pmesh = *gf.ParFESpace()->GetParMesh();
315  MPI_Comm comm = pmesh.GetComm();
316 
317  int num_procs, myid;
318  MPI_Comm_size(comm, &num_procs);
319  MPI_Comm_rank(comm, &myid);
320 
321  bool newly_opened = false;
322  int connection_failed;
323 
324  do
325  {
326  if (myid == 0)
327  {
328  if (!sock.is_open() || !sock)
329  {
330  sock.open(vishost, visport);
331  sock.precision(8);
332  newly_opened = true;
333  }
334  sock << "solution\n";
335  }
336 
337  pmesh.PrintAsOne(sock);
338  gf.SaveAsOne(sock);
339 
340  if (myid == 0 && newly_opened)
341  {
342  sock << "window_title '" << title << "'\n"
343  << "window_geometry "
344  << x << " " << y << " " << w << " " << h << "\n";
345  if ( keys ) { sock << "keys " << keys << "\n"; }
346  else { sock << "keys maaAc"; }
347  if ( vec ) { sock << "vvv"; }
348  sock << endl;
349  }
350 
351  if (myid == 0)
352  {
353  connection_failed = !sock && !newly_opened;
354  }
355  MPI_Bcast(&connection_failed, 1, MPI_INT, 0, comm);
356  }
357  while (connection_failed);
358 }
359 
360 } // namespace common
361 
362 } // namespace mfem
363 
364 #endif
const char vishost[]
void SetTol(double tol)
Definition: hypre.cpp:3995
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1032
IrrotationalProjector(ParFiniteElementSpace &H1FESpace, ParFiniteElementSpace &HCurlFESpace, const int &irOrder, ParBilinearForm *s0=NULL, ParMixedBilinearForm *weakDiv=NULL, ParDiscreteGradOperator *grad=NULL)
Definition: pfem_extras.cpp:98
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
virtual void SetIntRule(const IntegrationRule *ir)
Prescribe a fixed IntegrationRule to use (when ir != NULL) or let the integrator choose (when ir == N...
Definition: nonlininteg.hpp:43
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
const FiniteElementCollection * fec
Associated FE collection (not owned).
Definition: fespace.hpp:108
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:4015
void AddDomainInterpolator(DiscreteInterpolator *di)
Adds a domain interpolator. Assumes ownership of di.
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1579
ParDiscreteCurlOperator(ParFiniteElementSpace *dfes, ParFiniteElementSpace *rfes)
Definition: pfem_extras.cpp:83
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:81
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1659
void Assemble(int skip_zeros=1)
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void Assemble(int skip_zeros=1)
Assemble the local matrix.
const int visport
L2_ParFESpace(ParMesh *m, const int p, const int space_dim, int vdim=1, int order=Ordering::byNODES)
Definition: pfem_extras.cpp:61
void SetMaxIter(int max_iter)
Definition: hypre.cpp:4005
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:341
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:35
bool is_open()
True if the socketstream is open, false otherwise.
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
PCG solver in hypre.
Definition: hypre.hpp:1204
ParDiscreteGradOperator(ParFiniteElementSpace *dfes, ParFiniteElementSpace *rfes)
Definition: pfem_extras.cpp:76
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:899
ND_ParFESpace(ParMesh *m, const int p, const int space_dim, int vdim=1, int order=Ordering::byNODES)
Definition: pfem_extras.cpp:37
RT_ParFESpace(ParMesh *m, const int p, const int space_dim, int vdim=1, int order=Ordering::byNODES)
Definition: pfem_extras.cpp:49
Class for parallel bilinear form using different test and trial FE spaces.
void VisualizeMesh(socketstream &sock, const char *vishost, int visport, Mesh &mesh, const char *title, int x, int y, int w, int h, const char *keys)
Definition: fem_extras.cpp:59
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
ParDiscreteDivOperator(ParFiniteElementSpace *dfes, ParFiniteElementSpace *rfes)
Definition: pfem_extras.cpp:90
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:4020
Class for parallel bilinear form.
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:415
Vector data type.
Definition: vector.hpp:60
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)
Definition: fem_extras.cpp:92
void PrintAsOne(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4914
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:4043
Class for parallel meshes.
Definition: pmesh.hpp:32
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
DivergenceFreeProjector(ParFiniteElementSpace &H1FESpace, ParFiniteElementSpace &HCurlFESpace, const int &irOrder, ParBilinearForm *s0=NULL, ParMixedBilinearForm *weakDiv=NULL, ParDiscreteGradOperator *grad=NULL)
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:288
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113