MFEM  v3.3.2
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
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 miniapps
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 
127 
128  diffInteg->SetIntRule(ir);
129  wdivInteg->SetIntRule(ir);
130 
131  if ( s0 == NULL )
132  {
133  s0_ = new ParBilinearForm(H1FESpace_);
134  s0_->AddDomainIntegrator(diffInteg);
135  s0_->Assemble();
136  s0_->Finalize();
137  S0_ = new HypreParMatrix;
138  }
139  if ( weakDiv_ == NULL )
140  {
141  weakDiv_ = new ParMixedBilinearForm(HCurlFESpace_, H1FESpace_);
142  weakDiv_->AddDomainIntegrator(wdivInteg);
143  weakDiv_->Assemble();
144  weakDiv_->Finalize();
145  }
146  if ( grad_ == NULL )
147  {
148  grad_ = new ParDiscreteGradOperator(H1FESpace_, HCurlFESpace_);
149  grad_->Assemble();
150  grad_->Finalize();
151  }
152 
153  psi_ = new ParGridFunction(H1FESpace_);
154  xDiv_ = new ParGridFunction(H1FESpace_);
155 }
156 
158 {
159  delete psi_;
160  delete xDiv_;
161 
162  delete amg_;
163  delete pcg_;
164 
165  delete S0_;
166 
167  delete s0_;
168  delete weakDiv_;
169 }
170 
171 void
172 IrrotationalProjector::InitSolver() const
173 {
174  delete pcg_;
175  delete amg_;
176 
177  amg_ = new HypreBoomerAMG(*S0_);
178  amg_->SetPrintLevel(0);
179  pcg_ = new HyprePCG(*S0_);
180  pcg_->SetTol(1e-14);
181  pcg_->SetMaxIter(200);
182  pcg_->SetPrintLevel(0);
183  pcg_->SetPreconditioner(*amg_);
184 }
185 
186 void
188 {
189  // Compute the divergence of x
190  weakDiv_->Mult(x,*xDiv_); *xDiv_ *= -1.0;
191 
192  // Apply essential BC and form linear system
193  *psi_ = 0.0;
194  s0_->FormLinearSystem(ess_bdr_tdofs_, *psi_, *xDiv_, *S0_, Psi_, RHS_);
195 
196  // Solve the linear system for Psi
197  if ( pcg_ == NULL ) { this->InitSolver(); }
198  pcg_->Mult(RHS_, Psi_);
199 
200  // Compute the parallel grid function correspoinding to Psi
201  s0_->RecoverFEMSolution(Psi_, *xDiv_, *psi_);
202 
203  // Compute the irrotational portion of x
204  grad_->Mult(*psi_, y);
205 }
206 
207 void
209 {
210  delete pcg_; pcg_ = NULL;
211  delete amg_; amg_ = NULL;
212  delete S0_; S0_ = new HypreParMatrix;
213 
214  psi_->Update();
215  xDiv_->Update();
216 
217  if ( ownsS0_ )
218  {
219  s0_->Update();
220  s0_->Assemble();
221  s0_->Finalize();
222  }
223  if ( ownsWeakDiv_ )
224  {
225  weakDiv_->Update();
226  weakDiv_->Assemble();
227  weakDiv_->Finalize();
228  }
229  if ( ownsGrad_ )
230  {
231  grad_->Update();
232  grad_->Assemble();
233  grad_->Finalize();
234  }
235 
236  H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
237 }
238 
241  ParFiniteElementSpace & HCurlFESpace,
242  const int & irOrder,
243  ParBilinearForm * s0,
244  ParMixedBilinearForm * weakDiv,
246  : IrrotationalProjector(H1FESpace,HCurlFESpace, irOrder, s0, weakDiv, grad)
247 {}
248 
250 {}
251 
252 void
254 {
255  this->IrrotationalProjector::Mult(x, y);
256  y -= x;
257  y *= -1.0;
258 }
259 
260 void
262 {
264 }
265 
266 void VisualizeField(socketstream &sock, const char *vishost, int visport,
267  ParGridFunction &gf, const char *title,
268  int x, int y, int w, int h, bool vec)
269 {
270  ParMesh &pmesh = *gf.ParFESpace()->GetParMesh();
271  MPI_Comm comm = pmesh.GetComm();
272 
273  int num_procs, myid;
274  MPI_Comm_size(comm, &num_procs);
275  MPI_Comm_rank(comm, &myid);
276 
277  bool newly_opened = false;
278  int connection_failed;
279 
280  do
281  {
282  if (myid == 0)
283  {
284  if (!sock.is_open() || !sock)
285  {
286  sock.open(vishost, visport);
287  sock.precision(8);
288  newly_opened = true;
289  }
290  sock << "solution\n";
291  }
292 
293  pmesh.PrintAsOne(sock);
294  gf.SaveAsOne(sock);
295 
296  if (myid == 0 && newly_opened)
297  {
298  sock << "window_title '" << title << "'\n"
299  << "window_geometry "
300  << x << " " << y << " " << w << " " << h << "\n"
301  << "keys maaAc";
302  if ( vec ) { sock << "vvv"; }
303  sock << endl;
304  }
305 
306  if (myid == 0)
307  {
308  connection_failed = !sock && !newly_opened;
309  }
310  MPI_Bcast(&connection_failed, 1, MPI_INT, 0, comm);
311  }
312  while (connection_failed);
313 }
314 
315 } // namespace miniapps
316 
317 } // namespace mfem
318 
319 #endif
void SetTol(double tol)
Definition: hypre.cpp:2088
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:529
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
RT_ParFESpace(ParMesh *m, const int p, const int space_dim, int vdim=1, int order=Ordering::byNODES)
Definition: pfem_extras.cpp:49
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:861
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
L2_ParFESpace(ParMesh *m, const int p, const int space_dim, int vdim=1, int order=Ordering::byNODES)
Definition: pfem_extras.cpp:61
const Geometry::Type geom
void SaveAsOne(std::ostream &out=mfem::out)
Merge the local grid functions.
Definition: pgridfunc.cpp:376
Abstract parallel finite element space.
Definition: pfespace.hpp:31
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
const FiniteElementCollection * fec
Definition: fespace.hpp:66
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2103
ND_ParFESpace(ParMesh *m, const int p, const int space_dim, int vdim=1, int order=Ordering::byNODES)
Definition: pfem_extras.cpp:37
void AddDomainInterpolator(DiscreteInterpolator *di)
virtual void Update(FiniteElementSpace *nfes=NULL)
The BoomerAMG solver in hypre.
Definition: hypre.hpp:783
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:65
void PrintAsOne(std::ostream &out=mfem::out)
Definition: pmesh.cpp:3527
void SetPrintLevel(int print_level)
Definition: hypre.hpp:824
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, bool vec)
Definition: fem_extras.cpp:59
void Assemble(int skip_zeros=1)
void Assemble(int skip_zeros=1)
Assemble the local matrix.
ParDiscreteGradOperator(ParFiniteElementSpace *dfes, ParFiniteElementSpace *rfes)
Definition: pfem_extras.cpp:76
void SetMaxIter(int max_iter)
Definition: hypre.cpp:2093
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:239
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:22
MPI_Comm GetComm() const
Definition: pmesh.hpp:117
PCG solver in hypre.
Definition: hypre.hpp:643
DivergenceFreeProjector(ParFiniteElementSpace &H1FESpace, ParFiniteElementSpace &HCurlFESpace, const int &irOrder, ParBilinearForm *s0=NULL, ParMixedBilinearForm *weakDiv=NULL, ParDiscreteGradOperator *grad=NULL)
ParDiscreteDivOperator(ParFiniteElementSpace *dfes, ParFiniteElementSpace *rfes)
Definition: pfem_extras.cpp:90
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Class for parallel bilinear form using different test and trial FE spaces.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:2108
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:37
Class for parallel bilinear form.
int open(const char hostname[], int port)
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:297
Vector data type.
Definition: vector.hpp:41
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:146
IrrotationalProjector(ParFiniteElementSpace &H1FESpace, ParFiniteElementSpace &HCurlFESpace, const int &irOrder, ParBilinearForm *s0=NULL, ParMixedBilinearForm *weakDiv=NULL, ParDiscreteGradOperator *grad=NULL)
Definition: pfem_extras.cpp:98
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:175
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
virtual void Assemble(int skip_zeros=1)
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:2129
Class for parallel meshes.
Definition: pmesh.hpp:29
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:343
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:195
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:76