MFEM  v4.6.0
Finite element discretization library
extrapolator.cpp
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 #include "extrapolator.hpp"
13 #include "../common/mfem-common.hpp"
14 #include "marking.hpp"
15 
16 using namespace std;
17 
18 namespace mfem
19 {
20 
21 const char vishost[] = "localhost";
22 const int visport = 19916;
23 int wsize = 350; // glvis window size
24 
25 AdvectionOper::AdvectionOper(Array<bool> &zones, ParBilinearForm &Mbf,
26  ParBilinearForm &Kbf, const Vector &rhs)
27  : TimeDependentOperator(Mbf.Size()),
28  active_zones(zones),
29  M(Mbf), K(Kbf), K_mat(NULL), b(rhs),
30  lo_solver(NULL), lumpedM(NULL)
31 {
32  K_mat = K.ParallelAssemble(&K.SpMat());
33 
34  ParBilinearForm M_Lump(M.ParFESpace());
35  lumpedM = new Vector;
36  M_Lump.AddDomainIntegrator(new LumpedIntegrator(new MassIntegrator));
37  M_Lump.Assemble();
38  M_Lump.Finalize();
39  M_Lump.SpMat().GetDiag(*lumpedM);
40  lo_solver = new DiscreteUpwindLOSolver(*M.ParFESpace(),
41  K.SpMat(), *lumpedM);
42 }
43 
45 {
46  delete lo_solver;
47  delete lumpedM;
48  delete K_mat;
49 }
50 
51 void AdvectionOper::Mult(const Vector &x, Vector &dx) const
52 {
53  ParFiniteElementSpace &pfes = *M.ParFESpace();
54  const int NE = pfes.GetNE();
55  const int nd = pfes.GetFE(0)->GetDof();
56  Array<int> dofs(nd);
57 
58  if (adv_mode == LO)
59  {
60  lo_solver->CalcLOSolution(x, b, dx);
61  for (int k = 0; k < NE; k++)
62  {
63  pfes.GetElementDofs(k, dofs);
64  if (active_zones[k] == false)
65  {
66  dx.SetSubVector(dofs, 0.0);
67  continue;
68  }
69  }
70  return;
71  }
72 
73  MFEM_VERIFY(adv_mode == HO, "Wrong input for advection mode (-dg).");
74 
75  Vector rhs(x.Size());
76  K_mat->Mult(x, rhs);
77  rhs += b;
78 
79  DenseMatrix M_loc(nd);
80  DenseMatrixInverse M_loc_inv(&M_loc);
81  Vector rhs_loc(nd), dx_loc(nd);
82  for (int k = 0; k < NE; k++)
83  {
84  pfes.GetElementDofs(k, dofs);
85 
86  if (active_zones[k] == false)
87  {
88  dx.SetSubVector(dofs, 0.0);
89  continue;
90  }
91 
92  rhs.GetSubVector(dofs, rhs_loc);
93  M.SpMat().GetSubMatrix(dofs, dofs, M_loc);
94  M_loc_inv.Factor();
95  M_loc_inv.Mult(rhs_loc, dx_loc);
96  dx.SetSubVector(dofs, dx_loc);
97  }
98 }
99 
100 void AdvectionOper::ComputeElementsMinMax(const ParGridFunction &gf,
101  Vector &el_min, Vector &el_max) const
102 {
103  ParFiniteElementSpace &pfes = *gf.ParFESpace();
104  const int NE = pfes.GetNE(), ndof = pfes.GetFE(0)->GetDof();
105  for (int k = 0; k < NE; k++)
106  {
107  el_min(k) = numeric_limits<double>::infinity();
108  el_max(k) = -numeric_limits<double>::infinity();
109 
110  for (int i = 0; i < ndof; i++)
111  {
112  el_min(k) = min(el_min(k), gf(k*ndof + i));
113  el_max(k) = max(el_max(k), gf(k*ndof + i));
114  }
115  }
116 }
117 
118 void AdvectionOper::ComputeBounds(const ParFiniteElementSpace &pfes,
119  const Vector &el_min, const Vector &el_max,
120  Vector &dof_min, Vector &dof_max) const
121 {
122  ParMesh *pmesh = pfes.GetParMesh();
123  L2_FECollection fec_bounds(0, pmesh->Dimension());
124  ParFiniteElementSpace pfes_bounds(pmesh, &fec_bounds);
125  ParGridFunction el_min_gf(&pfes_bounds), el_max_gf(&pfes_bounds);
126  const int NE = pmesh->GetNE(), ndofs = dof_min.Size() / NE;
127 
128  el_min_gf = el_min;
129  el_max_gf = el_max;
130 
131  el_min_gf.ExchangeFaceNbrData(); el_max_gf.ExchangeFaceNbrData();
132  const Vector &min_nbr = el_min_gf.FaceNbrData();
133  const Vector &max_nbr = el_max_gf.FaceNbrData();
134  const Table &el_to_el = pmesh->ElementToElementTable();
135  Array<int> face_nbr_el;
136  for (int k = 0; k < NE; k++)
137  {
138  double k_min = el_min_gf(k), k_max = el_max_gf(k);
139 
140  el_to_el.GetRow(k, face_nbr_el);
141  for (int n = 0; n < face_nbr_el.Size(); n++)
142  {
143  if (face_nbr_el[n] < NE)
144  {
145  // Local neighbor.
146  k_min = std::min(k_min, el_min_gf(face_nbr_el[n]));
147  k_max = std::max(k_max, el_max_gf(face_nbr_el[n]));
148  }
149  else
150  {
151  // MPI face neighbor.
152  k_min = std::min(k_min, min_nbr(face_nbr_el[n] - NE));
153  k_max = std::max(k_max, max_nbr(face_nbr_el[n] - NE));
154  }
155  }
156 
157  for (int j = 0; j < ndofs; j++)
158  {
159  dof_min(k*ndofs + j) = k_min;
160  dof_max(k*ndofs + j) = k_max;
161  }
162  }
163 }
164 
166  const ParGridFunction &input,
167  const double time_period,
168  ParGridFunction &xtrap)
169 {
170  ParMesh &pmesh = *input.ParFESpace()->GetParMesh();
171  const int order = input.ParFESpace()->GetOrder(0),
172  dim = pmesh.Dimension(), NE = pmesh.GetNE();
173 
174  // Get a ParGridFunction and mark elements.
175  H1_FECollection fec(order, dim);
176  ParFiniteElementSpace pfes_H1(&pmesh, &fec);
177  ParGridFunction ls_gf(&pfes_H1);
178  ls_gf.ProjectCoefficient(level_set);
179  if (visualization)
180  {
181  socketstream sock1, sock2;
182  common::VisualizeField(sock1, vishost, visport, ls_gf,
183  "Domain level set", 0, 0, wsize, wsize,
184  "rRjlmm********A");
185  common::VisualizeField(sock2, vishost, visport, input,
186  "Input u", 0, wsize+60, wsize, wsize,
187  "rRjlmm********A");
188  MPI_Barrier(pmesh.GetComm());
189  }
190  // Mark elements.
191  Array<int> elem_marker;
192  ShiftedFaceMarker marker(pmesh, pfes_H1, false);
193  ls_gf.ExchangeFaceNbrData();
194  marker.MarkElements(ls_gf, elem_marker);
195 
196  // The active zones are where we extrapolate (where the PDE is solved).
197  Array<bool> active_zones(NE);
198  for (int k = 0; k < NE; k++)
199  {
200  // Extrapolation is done in zones that are CUT or OUTSIDE.
201  active_zones[k] =
202  (elem_marker[k] == ShiftedFaceMarker::INSIDE) ? false : true;
203  }
204 
205  // Setup a VectorCoefficient for n = - grad_ls / |grad_ls|.
206  // The sign makes it point out of the known region.
207  // The coefficient must be continuous to have well-defined transport.
208  LevelSetNormalGradCoeff ls_n_coeff_L2(ls_gf);
209  ParFiniteElementSpace pfes_H1_vec(&pmesh, &fec, dim);
210  ParGridFunction lsn_gf(&pfes_H1_vec);
211  ls_gf.ExchangeFaceNbrData();
212  lsn_gf.ProjectDiscCoefficient(ls_n_coeff_L2, GridFunction::ARITHMETIC);
213  VectorGridFunctionCoefficient ls_n_coeff(&lsn_gf);
214 
215  // Initial solution.
216  // Trim to the known values (only elements inside the known region).
217  Array<int> dofs;
218  L2_FECollection fec_L2(order, dim);
219  ParFiniteElementSpace pfes_L2(&pmesh, &fec_L2);
220  ParGridFunction u(&pfes_L2), vis_marking(&pfes_L2);
221  u.ProjectGridFunction(input);
222  for (int k = 0; k < NE; k++)
223  {
224  pfes_L2.GetElementDofs(k, dofs);
225  if (elem_marker[k] != ShiftedFaceMarker::INSIDE)
226  { u.SetSubVector(dofs, 0.0); }
227  vis_marking.SetSubVector(dofs, elem_marker[k]);
228  }
229  if (visualization)
230  {
231  socketstream sock1, sock2;
233  "Fixed (known) u values", wsize, 0,
234  wsize, wsize, "rRjlmm********A");
235  common::VisualizeField(sock2, vishost, visport, vis_marking,
236  "Element markings", 0, 2*wsize+60,
237  wsize, wsize, "rRjlmm********A");
238  }
239 
240  // Normal derivative function.
241  ParGridFunction n_grad_u(&pfes_L2);
242  NormalGradCoeff n_grad_u_coeff(u, ls_n_coeff);
243  n_grad_u.ProjectCoefficient(n_grad_u_coeff);
244  if (visualization && xtrap_degree >= 1)
245  {
246  socketstream sock;
247  common::VisualizeField(sock, vishost, visport, n_grad_u,
248  "n.grad(u)", 2*wsize, 0, wsize, wsize,
249  "rRjlmm********A");
250  }
251 
252  // 2nd normal derivative function.
253  ParGridFunction n_grad_n_grad_u(&pfes_L2);
254  NormalGradCoeff n_grad_n_grad_u_coeff(n_grad_u, ls_n_coeff);
255  n_grad_n_grad_u.ProjectCoefficient(n_grad_n_grad_u_coeff);
256  if (visualization && xtrap_degree == 2)
257  {
258  socketstream sock;
259  common::VisualizeField(sock, vishost, visport, n_grad_n_grad_u,
260  "n.grad(n.grad(u))", 3*wsize, 0, wsize, wsize,
261  "rRjmm********A");
262  }
263 
264  ParBilinearForm lhs_bf(&pfes_L2), rhs_bf(&pfes_L2);
265  lhs_bf.AddDomainIntegrator(new MassIntegrator);
266  const double alpha = -1.0;
267  rhs_bf.AddDomainIntegrator(new ConvectionIntegrator(ls_n_coeff, alpha));
268  auto trace_i = new NonconservativeDGTraceIntegrator(ls_n_coeff, alpha);
269  rhs_bf.AddInteriorFaceIntegrator(trace_i);
270  rhs_bf.KeepNbrBlock(true);
271 
272  ls_gf.ExchangeFaceNbrData();
273  lhs_bf.Assemble();
274  lhs_bf.Finalize();
275  rhs_bf.Assemble(0);
276  rhs_bf.Finalize(0);
277 
278  // Compute a CFL time step.
279  double h_min = std::numeric_limits<double>::infinity();
280  for (int k = 0; k < NE; k++)
281  {
282  h_min = std::min(h_min, pmesh.GetElementSize(k));
283  }
284  MPI_Allreduce(MPI_IN_PLACE, &h_min, 1, MPI_DOUBLE, MPI_MIN, pmesh.GetComm());
285  // The propagation speed is 1.
286  double dt = 0.25 * h_min / order / 1.0;
287  double half_dt = 0.5 * dt;
289  {
290  dt = half_dt;
291  }
292 
293  // Time loops.
294  Vector rhs(pfes_L2.GetVSize());
295  AdvectionOper adv_oper(active_zones, lhs_bf, rhs_bf, rhs);
296  adv_oper.adv_mode = advection_mode;
297  RK2Solver ode_solver(1.0);
298  ode_solver.Init(adv_oper);
299 
300  if (xtrap_degree == 0)
301  {
302  // Constant extrapolation of u (always LO).
303  rhs = 0.0;
304  adv_oper.adv_mode = AdvectionOper::LO;
305  TimeLoop(u, ode_solver, time_period, half_dt,
306  wsize, "Extrap const u -- LO");
307  xtrap.ProjectGridFunction(u);
308  return;
309  }
310 
311  std::string mode_text = "HO";
312  if (advection_mode == AdvectionOper::LO) { mode_text = "LO"; }
313 
314  MFEM_VERIFY(xtrap_degree == 1 || xtrap_degree == 2, "Wrong order input.");
315  if (xtrap_type == ASLAM)
316  {
317  if (xtrap_degree == 1)
318  {
319  // Constant extrapolation of [n.grad_u] (always LO).
320  rhs = 0.0;
321  adv_oper.adv_mode = AdvectionOper::LO;
322  TimeLoop(n_grad_u, ode_solver, time_period, half_dt,
323  2*wsize, "Extrap const n.grad(u) -- Aslam -- LO");
324 
325  adv_oper.adv_mode = advection_mode;
326 
327  // Linear extrapolation of u.
328  lhs_bf.Mult(n_grad_u, rhs);
329  TimeLoop(u, ode_solver, time_period, dt,
330  wsize, "Extrap linear u -- Aslam -- " + mode_text);
331  }
332 
333  if (xtrap_degree == 2)
334  {
335  // Constant extrapolation of [n.grad(n.grad(u))] (always LO).
336  rhs = 0.0;
337  adv_oper.adv_mode = AdvectionOper::LO;
338  TimeLoop(n_grad_n_grad_u, ode_solver, time_period, half_dt,
339  3*wsize, "Extrap const n.grad(n.grad(u)) -- Aslam -- LO");
340 
341  adv_oper.adv_mode = advection_mode;
342 
343  // Linear extrapolation of [n.grad_u].
344  lhs_bf.Mult(n_grad_n_grad_u, rhs);
345  TimeLoop(n_grad_u, ode_solver, time_period, dt,
346  2*wsize, "Extrap linear n.grad(u) -- Aslam -- " + mode_text);
347 
348  // Quadratic extrapolation of u.
349  lhs_bf.Mult(n_grad_u, rhs);
350  TimeLoop(u, ode_solver, time_period, dt,
351  wsize, "Extrap quadratic u -- Aslam -- " + mode_text);
352  }
353  }
354  else if (xtrap_type == BOCHKOV)
355  {
356  if (xtrap_degree == 1)
357  {
358  // Constant extrapolation of all grad(u) components (always LO).
359  rhs = 0.0;
360  adv_oper.adv_mode = AdvectionOper::LO;
361  ParGridFunction grad_u_0(&pfes_L2), grad_u_1(&pfes_L2);
362  GradComponentCoeff grad_u_0_coeff(u, 0), grad_u_1_coeff(u, 1);
363  grad_u_0.ProjectCoefficient(grad_u_0_coeff);
364  grad_u_1.ProjectCoefficient(grad_u_1_coeff);
365  TimeLoop(grad_u_0, ode_solver, time_period, half_dt,
366  2*wsize, "Extrap const du_dx -- Bochkov -- LO");
367  TimeLoop(grad_u_1, ode_solver, time_period, half_dt,
368  3*wsize, "Extrap const du_dy -- Bochkov -- LO");
369 
370  adv_oper.adv_mode = advection_mode;
371 
372  // Linear extrapolation of u.
373  ParLinearForm rhs_lf(&pfes_L2);
374  NormalGradComponentCoeff grad_u_n(grad_u_0, grad_u_1, ls_n_coeff);
375  rhs_lf.AddDomainIntegrator(new DomainLFIntegrator(grad_u_n));
376  rhs_lf.Assemble();
377  rhs = rhs_lf;
378  TimeLoop(u, ode_solver, time_period, dt,
379  wsize, "Extrap linear u -- Bochkov -- " + mode_text);
380  }
381 
382  if (xtrap_degree == 2)
383  {
384  MFEM_ABORT("Quadratic Bochkov method is not implemented.");
385  }
386  }
387  else { MFEM_ABORT("Wrong input for extrapolation type (-et)."); }
388 
389  xtrap.ProjectGridFunction(u);
390 }
391 
392 // Errors in cut elements.
394  const ParGridFunction &exact,
395  const ParGridFunction &xtrap,
396  double &err_L1, double &err_L2,
397  double &err_LI)
398 {
399  ParMesh &pmesh = *exact.ParFESpace()->GetParMesh();
400  const int order = exact.ParFESpace()->GetOrder(0),
401  dim = pmesh.Dimension(), NE = pmesh.GetNE();
402 
403  // Get a ParGridFunction and mark elements.
404  H1_FECollection fec(order, dim);
405  ParFiniteElementSpace pfes_H1(&pmesh, &fec);
406  ParGridFunction ls_gf(&pfes_H1);
407  ls_gf.ProjectCoefficient(level_set);
408  // Mark elements.
409  Array<int> elem_marker;
410  ShiftedFaceMarker marker(pmesh, pfes_H1, false);
411  ls_gf.ExchangeFaceNbrData();
412  marker.MarkElements(ls_gf, elem_marker);
413 
414  Vector errors_L1(NE), errors_L2(NE), errors_LI(NE);
415  GridFunctionCoefficient exact_coeff(&exact);
416 
417  xtrap.ComputeElementL1Errors(exact_coeff, errors_L1);
418  xtrap.ComputeElementL2Errors(exact_coeff, errors_L2);
419  xtrap.ComputeElementMaxErrors(exact_coeff, errors_LI);
420  err_L1 = 0.0, err_L2 = 0.0, err_LI = 0.0;
421  double cut_volume = 0.0;
422  for (int k = 0; k < NE; k++)
423  {
424  if (elem_marker[k] == ShiftedFaceMarker::CUT)
425  {
426  err_L1 += errors_L1(k);
427  err_L2 += errors_L2(k);
428  err_LI = std::max(err_LI, errors_LI(k));
429  cut_volume += pmesh.GetElementVolume(k);
430  }
431  }
432  MPI_Comm comm = pmesh.GetComm();
433  MPI_Allreduce(MPI_IN_PLACE, &err_L1, 1, MPI_DOUBLE, MPI_SUM, comm);
434  MPI_Allreduce(MPI_IN_PLACE, &err_L2, 1, MPI_DOUBLE, MPI_SUM, comm);
435  MPI_Allreduce(MPI_IN_PLACE, &err_LI, 1, MPI_DOUBLE, MPI_MAX, comm);
436  MPI_Allreduce(MPI_IN_PLACE, &cut_volume, 1, MPI_DOUBLE, MPI_SUM, comm);
437  err_L1 /= cut_volume;
438  err_L2 /= cut_volume;
439 }
440 
441 void Extrapolator::TimeLoop(ParGridFunction &sltn, ODESolver &ode_solver,
442  double t_final, double dt,
443  int vis_x_pos, std::string vis_name)
444 {
445  socketstream sock;
446 
447  const int myid = sltn.ParFESpace()->GetMyRank();
448  bool done = false;
449  double t = 0.0;
450  for (int ti = 0; !done;)
451  {
452  double dt_real = min(dt, t_final - t);
453  ode_solver.Step(sltn, t, dt_real);
454  ti++;
455 
456  done = (t >= t_final - 1e-8*dt);
457  if (done || ti % vis_steps == 0)
458  {
459  if (myid == 0)
460  {
461  cout << vis_name+" / time step: " << ti << ", time: " << t << endl;
462  }
463  if (visualization)
464  {
466  vis_name.c_str(), vis_x_pos, wsize+60,
467  wsize, wsize, "rRjlmm********A");
468  MPI_Barrier(sltn.ParFESpace()->GetComm());
469  }
470  }
471  }
472 }
473 
474 
475 
477  const SparseMatrix &adv,
478  const Vector &Mlump)
479  : pfes(space), K(adv), D(adv), K_smap(), M_lumped(Mlump)
480 {
481  // Assuming it is finalized.
482  const int *I = K.GetI(), *J = K.GetJ(), n = K.Size();
483  K_smap.SetSize(I[n]);
484  for (int row = 0, j = 0; row < n; row++)
485  {
486  for (int end = I[row+1]; j < end; j++)
487  {
488  int col = J[j];
489  // Find the offset, _j, of the (col,row) entry and store it in smap[j].
490  for (int _j = I[col], _end = I[col+1]; true; _j++)
491  {
492  MFEM_VERIFY(_j != _end, "Can't find the symmetric entry!");
493 
494  if (J[_j] == row) { K_smap[j] = _j; break; }
495  }
496  }
497  }
498 
500 }
501 
503  Vector &du) const
504 {
505  ParGridFunction u_gf(&pfes);
506  u_gf = u;
507  ApplyDiscreteUpwindMatrix(u_gf, du);
508 
509  const int s = du.Size();
510  for (int i = 0; i < s; i++)
511  {
512  du(i) = (du(i) + rhs(i)) / M_lumped(i);
513  }
514 }
515 
517 {
518  const int *I = K.HostReadI(), *J = K.HostReadJ(), n = K.Size();
519 
520  const double *K_data = K.HostReadData();
521 
522  double *D_data = D.HostReadWriteData();
524 
525  for (int i = 0, k = 0; i < n; i++)
526  {
527  double rowsum = 0.;
528  for (int end = I[i+1]; k < end; k++)
529  {
530  int j = J[k];
531  double kij = K_data[k];
532  double kji = K_data[K_smap[k]];
533  double dij = fmax(fmax(0.0,-kij),-kji);
534  D_data[k] = kij + dij;
535  D_data[K_smap[k]] = kji + dij;
536  if (i != j) { rowsum += dij; }
537  }
538  D(i,i) = K(i,i) - rowsum;
539  }
540 }
541 
543  Vector &du) const
544 {
545  const int s = u.Size();
546  const int *I = D.HostReadI(), *J = D.HostReadJ();
547  const double *D_data = D.HostReadData();
548 
549  u.ExchangeFaceNbrData();
550  const Vector &u_np = u.FaceNbrData();
551 
552  for (int i = 0; i < s; i++)
553  {
554  du(i) = 0.0;
555  for (int k = I[i]; k < I[i + 1]; k++)
556  {
557  int j = J[k];
558  double u_j = (j < s) ? u(j) : u_np[j - s];
559  double d_ij = D_data[k];
560  du(i) += d_ij * u_j;
561  }
562  }
563 }
564 
565 }
const char vishost[]
const double * HostReadData() const
Definition: sparsemat.hpp:277
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
enum mfem::Extrapolator::XtrapType xtrap_type
ParFiniteElementSpace * ParFESpace() const
Return the parallel FE space associated with the ParBilinearForm.
int Size() const
For backward compatibility, define Size() to be synonym of Height().
Definition: sparsemat.hpp:206
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
void KeepNbrBlock(bool knb=true)
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3943
Base abstract class for first order time dependent operators.
Definition: operator.hpp:316
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:227
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:23
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
Project a discontinuous vector coefficient as a grid function on a continuous finite element space...
Definition: pgridfunc.cpp:582
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3897
STL namespace.
virtual void ComputeElementMaxErrors(Coefficient &exsol, Vector &error, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:630
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
Definition: sparsemat.cpp:2886
DiscreteUpwindLOSolver(ParFiniteElementSpace &space, const SparseMatrix &adv, const Vector &Mlump)
double * HostReadWriteData()
Definition: sparsemat.hpp:281
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
const SparseMatrix & K
void Extrapolate(Coefficient &level_set, const ParGridFunction &input, const double time_period, ParGridFunction &xtrap)
Class for parallel linear form.
Definition: plinearform.hpp:26
double GetElementSize(ElementTransformation *T, int type=0)
Definition: mesh.cpp:83
Data type sparse matrix.
Definition: sparsemat.hpp:50
void CalcLOSolution(const Vector &u, const Vector &rhs, Vector &du) const
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
double b
Definition: lissajous.cpp:42
const int * HostReadJ() const
Definition: sparsemat.hpp:261
ParMesh * GetParMesh() const
Definition: pfespace.hpp:273
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void ComputeLocalErrors(Coefficient &level_set, const ParGridFunction &exact, const ParGridFunction &xtrap, double &err_L1, double &err_L2, double &err_LI)
const int visport
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:736
void ComputeDiscreteUpwindMatrix() const
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
virtual DofTransformation * GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: pfespace.cpp:470
int wsize
int wsize
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:534
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
void ProjectGridFunction(const GridFunction &src)
Project the src GridFunction to this GridFunction, both of which must be on the same mesh...
Definition: gridfunc.cpp:1849
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
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
string space
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
int * HostReadWriteJ()
Definition: sparsemat.hpp:265
virtual void ComputeElementL2Errors(Coefficient &exsol, Vector &error, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:624
const int * HostReadI() const
Definition: sparsemat.hpp:245
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:39
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
virtual void ComputeElementL1Errors(Coefficient &exsol, Vector &error, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:618
virtual void Mult(const Vector &x, Vector &dx) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
int dim
Definition: ex24.cpp:53
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
enum mfem::AdvectionOper::AdvectionMode adv_mode
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:44
Class for parallel bilinear form.
RefCoord t[3]
void ApplyDiscreteUpwindMatrix(ParGridFunction &u, Vector &du) const
MPI_Comm GetComm() const
Definition: pfespace.hpp:269
const double alpha
Definition: ex15.cpp:369
Vector data type.
Definition: vector.hpp:58
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
Vector coefficient defined by a vector GridFunction.
ParFiniteElementSpace & pfes
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
int * HostReadWriteI()
Definition: sparsemat.hpp:249
RefCoord s[3]
AdvectionOper::AdvectionMode advection_mode
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Class for parallel meshes.
Definition: pmesh.hpp:32
void MarkElements(const ParGridFunction &ls_func, Array< int > &elem_marker)
Definition: marking.cpp:17
double GetElementVolume(int i)
Definition: mesh.cpp:119
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1736
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
alpha (q . grad u, v)