MFEM  v3.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
pgridfunc.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.googlecode.com.
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 "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "fem.hpp"
17 #include <iostream>
18 #include <limits>
19 using namespace std;
20 
21 namespace mfem
22 {
23 
24 ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf, GridFunction *gf)
25 {
26  fes = pfes = pf;
27  SetDataAndSize(gf->GetData(), gf->Size());
28 }
29 
30 ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf, HypreParVector *tv)
31  : GridFunction(pf), pfes(pf)
32 {
33  Distribute(tv);
34 }
35 
36 ParGridFunction::ParGridFunction(ParMesh *pmesh, GridFunction *gf, int * partitioning)
37 {
38  // duplicate the FiniteElementCollection from 'gf'
40  fes = pfes = new ParFiniteElementSpace(pmesh, fec, gf->FESpace()->GetVDim(),
41  gf->FESpace()->GetOrdering());
42  SetSize(pfes->GetVSize());
43 
44  if(partitioning)
45  {
46  Array<int> gvdofs, lvdofs;
47  Vector lnodes;
48  int element_counter = 0;
49  Mesh & mesh(*gf->FESpace()->GetMesh());
50  int MyRank;
51  MPI_Comm_rank(pfes->GetComm(), &MyRank);
52  for (int i = 0; i < mesh.GetNE(); i++)
53  if (partitioning[i] == MyRank)
54  {
55  pfes->GetElementVDofs(element_counter, lvdofs);
56  gf->FESpace()->GetElementVDofs(i, gvdofs);
57  gf->GetSubVector(gvdofs, lnodes);
58  SetSubVector(lvdofs, lnodes);
59  element_counter++;
60  }
61  }
62 }
63 
65 {
68  pfes = f;
69 }
70 
72 {
74  GridFunction::Update(f, v, v_offset);
75  pfes = f;
76 }
77 
79 {
80  pfes->Dof_TrueDof_Matrix()->Mult(*tv, *this);
81 }
82 
84 {
85 #if 0
86  for (int i = 0; i < size; i++)
87  {
88  int tdof = pfes->GetLocalTDofNumber(i);
89  if (tdof >= 0)
90  tv(tdof) = (*this)(i);
91  }
92 #else
93  hypre_ParCSRMatrix *P = *pfes->Dof_TrueDof_Matrix();
94  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(P);
95  int *I = hypre_CSRMatrixI(diag) + 1;
96  int *J = hypre_CSRMatrixJ(diag);
97  for (int i = 0, j = 0; i < size; i++)
98  if (j < I[i])
99  tv(J[j++]) = (*this)(i);
100 #endif
101 }
102 
104 {
106  GetTrueDofs(*tv);
107  return tv;
108 }
109 
111 {
112  pfes->Dof_TrueDof_Matrix()->MultTranspose(*this, tv);
113  pfes->DivideByGroupSize(tv);
114 }
115 
117 {
118  pfes->Dof_TrueDof_Matrix()->MultTranspose(*this, tv);
119  pfes->DivideByGroupSize(tv);
120 }
121 
123 {
125  ParallelAverage(*tv);
126  return tv;
127 }
128 
130 {
131  pfes->Dof_TrueDof_Matrix()->MultTranspose(*this, tv);
132 }
133 
135 {
136  pfes->Dof_TrueDof_Matrix()->MultTranspose(*this, tv);
137 }
138 
140 {
142  ParallelAssemble(*tv);
143  return tv;
144 }
145 
147 {
149 
150  if (pfes->GetFaceNbrVSize() <= 0)
151  return;
152 
153  ParMesh *pmesh = pfes->GetParMesh();
154 
157 
158  int *send_offset = pfes->send_face_nbr_ldof.GetI();
159  int *send_ldof = pfes->send_face_nbr_ldof.GetJ();
160  int *recv_offset = pfes->face_nbr_gdof.GetI();
161  MPI_Comm MyComm = pfes->GetComm();
162 
163  int num_face_nbrs = pmesh->GetNFaceNeighbors();
164  MPI_Request *requests = new MPI_Request[2*num_face_nbrs];
165  MPI_Request *send_requests = requests;
166  MPI_Request *recv_requests = requests + num_face_nbrs;
167  MPI_Status *statuses = new MPI_Status[num_face_nbrs];
168 
169  for (int i = 0; i < send_data.Size(); i++)
170  send_data[i] = data[send_ldof[i]];
171 
172  for (int fn = 0; fn < num_face_nbrs; fn++)
173  {
174  int nbr_rank = pmesh->GetFaceNbrRank(fn);
175  int tag = 0;
176 
177  MPI_Isend(&send_data(send_offset[fn]),
178  send_offset[fn+1] - send_offset[fn],
179  MPI_DOUBLE, nbr_rank, tag, MyComm, &send_requests[fn]);
180 
181  MPI_Irecv(&face_nbr_data(recv_offset[fn]),
182  recv_offset[fn+1] - recv_offset[fn],
183  MPI_DOUBLE, nbr_rank, tag, MyComm, &recv_requests[fn]);
184  }
185 
186  MPI_Waitall(num_face_nbrs, send_requests, statuses);
187  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
188 
189  delete [] statuses;
190  delete [] requests;
191 }
192 
193 double ParGridFunction::GetValue(int i, const IntegrationPoint &ip, int vdim)
194  const
195 {
196  Array<int> dofs;
197  Vector DofVal, LocVec;
198  int nbr_el_no = i - pfes->GetParMesh()->GetNE();
199  if (nbr_el_no >= 0)
200  {
201  int fes_vdim = pfes->GetVDim();
202  pfes->GetFaceNbrElementVDofs(nbr_el_no, dofs);
203  if (fes_vdim > 1)
204  {
205  int s = dofs.Size()/fes_vdim;
206  Array<int> _dofs(&dofs[(vdim-1)*s], s);
207  face_nbr_data.GetSubVector(_dofs, LocVec);
208  DofVal.SetSize(s);
209  }
210  else
211  {
212  face_nbr_data.GetSubVector(dofs, LocVec);
213  DofVal.SetSize(dofs.Size());
214  }
215  pfes->GetFaceNbrFE(nbr_el_no)->CalcShape(ip, DofVal);
216  }
217  else
218  {
219  fes->GetElementDofs(i, dofs);
220  fes->DofsToVDofs(vdim-1, dofs);
221  DofVal.SetSize(dofs.Size());
222  fes->GetFE(i)->CalcShape(ip, DofVal);
223  GetSubVector(dofs, LocVec);
224  }
225 
226  return (DofVal * LocVec);
227 }
228 
230 {
231  DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
232 
233  if (delta_c == NULL)
234  {
236  }
237  else
238  {
239  double loc_integral, glob_integral;
240 
241  ProjectDeltaCoefficient(*delta_c, loc_integral);
242 
243  MPI_Allreduce(&loc_integral, &glob_integral, 1, MPI_DOUBLE, MPI_SUM,
244  pfes->GetComm());
245 
246  (*this) *= (delta_c->Scale() / glob_integral);
247  }
248 }
249 
250 void ParGridFunction::Save(std::ostream &out) const
251 {
252  for (int i = 0; i < size; i++)
253  if (pfes->GetDofSign(i) < 0)
254  data[i] = -data[i];
255 
256  GridFunction::Save(out);
257 
258  for (int i = 0; i < size; i++)
259  if (pfes->GetDofSign(i) < 0)
260  data[i] = -data[i];
261 }
262 
263 void ParGridFunction::SaveAsOne(std::ostream &out)
264 {
265  int i, p;
266 
267  MPI_Comm MyComm;
268  MPI_Status status;
269  int MyRank, NRanks;
270 
271  MyComm = pfes -> GetComm();
272 
273  MPI_Comm_size(MyComm, &NRanks);
274  MPI_Comm_rank(MyComm, &MyRank);
275 
276  double **values = new double*[NRanks];
277  int *nv = new int[NRanks];
278  int *nvdofs = new int[NRanks];
279  int *nedofs = new int[NRanks];
280  int *nfdofs = new int[NRanks];
281  int *nrdofs = new int[NRanks];
282 
283  values[0] = data;
284  nv[0] = pfes -> GetVSize();
285  nvdofs[0] = pfes -> GetNVDofs();
286  nedofs[0] = pfes -> GetNEDofs();
287  nfdofs[0] = pfes -> GetNFDofs();
288 
289  if (MyRank == 0)
290  {
291  pfes -> Save(out);
292  out << '\n';
293 
294  for (p = 1; p < NRanks; p++)
295  {
296  MPI_Recv(&nv[p], 1, MPI_INT, p, 455, MyComm, &status);
297  MPI_Recv(&nvdofs[p], 1, MPI_INT, p, 456, MyComm, &status);
298  MPI_Recv(&nedofs[p], 1, MPI_INT, p, 457, MyComm, &status);
299  MPI_Recv(&nfdofs[p], 1, MPI_INT, p, 458, MyComm, &status);
300  values[p] = new double[nv[p]];
301  MPI_Recv(values[p], nv[p], MPI_DOUBLE, p, 460, MyComm, &status);
302  }
303 
304  int vdim = pfes -> GetVDim();
305 
306  for (p = 0; p < NRanks; p++)
307  nrdofs[p] = nv[p]/vdim - nvdofs[p] - nedofs[p] - nfdofs[p];
308 
310  {
311  for (int d = 0; d < vdim; d++)
312  {
313  for (p = 0; p < NRanks; p++)
314  for (i = 0; i < nvdofs[p]; i++)
315  out << *values[p]++ << '\n';
316 
317  for (p = 0; p < NRanks; p++)
318  for (i = 0; i < nedofs[p]; i++)
319  out << *values[p]++ << '\n';
320 
321  for (p = 0; p < NRanks; p++)
322  for (i = 0; i < nfdofs[p]; i++)
323  out << *values[p]++ << '\n';
324 
325  for (p = 0; p < NRanks; p++)
326  for (i = 0; i < nrdofs[p]; i++)
327  out << *values[p]++ << '\n';
328  }
329  }
330  else
331  {
332  for (p = 0; p < NRanks; p++)
333  for (i = 0; i < nvdofs[p]; i++)
334  for (int d = 0; d < vdim; d++)
335  out << *values[p]++ << '\n';
336 
337  for (p = 0; p < NRanks; p++)
338  for (i = 0; i < nedofs[p]; i++)
339  for (int d = 0; d < vdim; d++)
340  out << *values[p]++ << '\n';
341 
342  for (p = 0; p < NRanks; p++)
343  for (i = 0; i < nfdofs[p]; i++)
344  for (int d = 0; d < vdim; d++)
345  out << *values[p]++ << '\n';
346 
347  for (p = 0; p < NRanks; p++)
348  for (i = 0; i < nrdofs[p]; i++)
349  for (int d = 0; d < vdim; d++)
350  out << *values[p]++ << '\n';
351  }
352 
353  for (p = 1; p < NRanks; p++)
354  {
355  values[p] -= nv[p];
356  delete [] values[p];
357  }
358  }
359  else
360  {
361  MPI_Send(&nv[0], 1, MPI_INT, 0, 455, MyComm);
362  MPI_Send(&nvdofs[0], 1, MPI_INT, 0, 456, MyComm);
363  MPI_Send(&nedofs[0], 1, MPI_INT, 0, 457, MyComm);
364  MPI_Send(&nfdofs[0], 1, MPI_INT, 0, 458, MyComm);
365  MPI_Send(data, nv[0], MPI_DOUBLE, 0, 460, MyComm);
366  }
367 
368  delete [] values;
369  delete [] nv;
370  delete [] nvdofs;
371  delete [] nedofs;
372  delete [] nfdofs;
373  delete [] nrdofs;
374 }
375 
376 double GlobalLpNorm(const double p, double loc_norm, MPI_Comm comm)
377 {
378  double glob_norm;
379 
380  if (p < numeric_limits<double>::infinity())
381  {
382  // negative quadrature weights may cause the error to be negative
383  if (loc_norm < 0.0)
384  loc_norm = -pow(-loc_norm, p);
385  else
386  loc_norm = pow(loc_norm, p);
387 
388  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
389 
390  if (glob_norm < 0.0)
391  glob_norm = -pow(-glob_norm, 1.0/p);
392  else
393  glob_norm = pow(glob_norm, 1.0/p);
394  }
395  else
396  {
397  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
398  }
399 
400  return glob_norm;
401 }
402 
403 }
404 
405 #endif
int GetNFaceNeighbors() const
Definition: pmesh.hpp:99
int Size() const
Logical size of the array.
Definition: array.hpp:108
int GetVSize() const
Definition: fespace.hpp:151
int * GetJ()
Definition: table.hpp:88
HypreParVector * NewTrueDofVector()
Definition: pfespace.hpp:127
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:26
void SaveAsOne(std::ostream &out=std::cout)
Merge the local grid functions.
Definition: pgridfunc.cpp:263
HypreParVector * ParallelAssemble() const
Returns a new vector assembled on the true dofs.
Definition: pgridfunc.cpp:139
void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff, double &integral)
Definition: gridfunc.cpp:906
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:248
ParFiniteElementSpace * pfes
Definition: pgridfunc.hpp:34
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:403
Delta function coefficient.
int Size() const
Returns the size of the vector.
Definition: vector.hpp:76
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:396
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: pgridfunc.cpp:193
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:250
void DivideByGroupSize(double *vec)
Scale a vector of true dofs.
Definition: pfespace.cpp:401
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:229
double * GetData() const
Definition: vector.hpp:80
int Size_of_connections() const
Definition: table.hpp:72
int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:572
void GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:689
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:1797
const FiniteElement * GetFaceNbrFE(int i) const
Definition: pfespace.cpp:695
int MultTranspose(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A^t * x + beta * y.
Definition: hypre.cpp:631
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:132
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:739
int GetLocalTDofNumber(int ldof)
Definition: pfespace.cpp:440
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:143
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:246
HypreParMatrix * Dof_TrueDof_Matrix()
The true dof-to-dof interpolation matrix.
Definition: pfespace.cpp:303
int GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:160
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:38
double GlobalLpNorm(const double p, double loc_norm, MPI_Comm comm)
Compute a global Lp norm from the local Lp norms computed by each processor.
Definition: pgridfunc.cpp:376
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
FiniteElementSpace * fes
FE space on which grid function lives.
Definition: gridfunc.hpp:30
void SetSubVector(const Array< int > &dofs, const Vector &elemvect)
Definition: vector.cpp:427
void GetElementVDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:117
FiniteElementCollection * fec
Used when the grid function is read from a file.
Definition: gridfunc.hpp:33
void Distribute(const Vector *tv)
Definition: pgridfunc.cpp:78
static FiniteElementCollection * New(const char *name)
Definition: fe_coll.cpp:38
virtual const char * Name() const
Definition: fe_coll.hpp:36
Class for integration point with weight.
Definition: intrules.hpp:25
HypreParVector * ParallelAverage() const
Returns a new vector averaged on the true dofs.
Definition: pgridfunc.cpp:122
void Destroy()
Destroy a vector.
Definition: vector.hpp:263
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:814
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:103
void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:967
int GetFaceNbrVSize() const
Definition: pfespace.hpp:160
void DofsToVDofs(Array< int > &dofs) const
Definition: fespace.cpp:29
Vector data type.
Definition: vector.hpp:29
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:162
int * GetI()
Definition: table.hpp:87
double * data
Definition: vector.hpp:34
Class for parallel meshes.
Definition: pmesh.hpp:27
int GetFaceNbrRank(int fn) const
Definition: pmesh.cpp:1082