MFEM  v4.5.2
Finite element discretization library
dist_solver.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 "dist_solver.hpp"
13 
14 namespace mfem
15 {
16 
17 void DiffuseField(ParGridFunction &field, int smooth_steps)
18 {
19  // Setup the Laplacian operator.
20  ParBilinearForm *Lap = new ParBilinearForm(field.ParFESpace());
22  Lap->Assemble();
23  Lap->Finalize();
24  HypreParMatrix *A = Lap->ParallelAssemble();
25 
26  HypreSmoother *S = new HypreSmoother(*A,0,smooth_steps);
27  S->iterative_mode = true;
28 
29  Vector tmp(A->Width());
30  field.SetTrueVector();
31  Vector fieldtrue = field.GetTrueVector();
32  tmp = 0.0;
33  S->Mult(tmp, fieldtrue);
34 
35  field.SetFromTrueDofs(fieldtrue);
36 
37  delete A;
38  delete S;
39  delete Lap;
40 }
41 
42 double AvgElementSize(ParMesh &pmesh)
43 {
44  // Compute average mesh size (assumes similar cells).
45  double dx, loc_area = 0.0;
46  for (int i = 0; i < pmesh.GetNE(); i++)
47  {
48  loc_area += pmesh.GetElementVolume(i);
49  }
50  double glob_area;
51  MPI_Allreduce(&loc_area, &glob_area, 1, MPI_DOUBLE,
52  MPI_SUM, pmesh.GetComm());
53  const int glob_zones = pmesh.GetGlobalNE();
54  switch (pmesh.GetElementBaseGeometry(0))
55  {
56  case Geometry::SEGMENT:
57  dx = glob_area / glob_zones; break;
58  case Geometry::SQUARE:
59  dx = sqrt(glob_area / glob_zones); break;
60  case Geometry::TRIANGLE:
61  dx = sqrt(2.0 * glob_area / glob_zones); break;
62  case Geometry::CUBE:
63  dx = pow(glob_area / glob_zones, 1.0/3.0); break;
65  dx = pow(6.0 * glob_area / glob_zones, 1.0/3.0); break;
66  default: MFEM_ABORT("Unknown zone type!"); dx = 0.0;
67  }
68 
69  return dx;
70 }
71 
73  ParGridFunction &dist_v)
74 {
75  ParFiniteElementSpace &pfes = *dist_s.ParFESpace();
76  MFEM_VERIFY(pfes.GetOrdering()==Ordering::byNODES,
77  "Only Ordering::byNODES is supported.");
78 
79  const int dim = pfes.GetMesh()->Dimension();
80  const int size = dist_s.Size();
81 
82  ParGridFunction der(&pfes);
83  Vector magn(size);
84  magn = 0.0;
85  for (int d = 0; d < dim; d++)
86  {
87  dist_s.GetDerivative(1, d, der);
88  for (int i = 0; i < size; i++)
89  {
90  magn(i) += der(i) * der(i);
91  // The vector must point towards the level zero set.
92  dist_v(i + d*size) = (dist_s(i) > 0.0) ? -der(i) : der(i);
93  }
94  }
95 
96  for (int i = 0; i < size; i++)
97  {
98  const double vec_magn = std::sqrt(magn(i) + 1e-12);
99  for (int d = 0; d < dim; d++)
100  {
101  dist_v(i + d*size) *= fabs(dist_s(i)) / vec_magn;
102  }
103  }
104 }
105 
107  ParGridFunction &distance)
108 {
109  ParFiniteElementSpace &pfes = *distance.ParFESpace();
110  MFEM_VERIFY(pfes.GetVDim() == pfes.GetMesh()->Dimension(),
111  "This function expects a vector ParGridFunction!");
112 
113  ParFiniteElementSpace pfes_s(pfes.GetParMesh(), pfes.FEColl());
114  ParGridFunction dist_s(&pfes_s);
115  ComputeScalarDistance(zero_level_set, dist_s);
116  ScalarDistToVector(dist_s, distance);
117 }
118 
119 
121  ParGridFunction &distance)
122 {
123  ParFiniteElementSpace &pfes = *distance.ParFESpace();
124 
125  auto check_h1 = dynamic_cast<const H1_FECollection *>(pfes.FEColl());
126  MFEM_VERIFY(check_h1 && pfes.GetVDim() == 1,
127  "This solver supports only scalar H1 spaces.");
128 
129  // Compute average mesh size (assumes similar cells).
130  ParMesh &pmesh = *pfes.GetParMesh();
131 
132  // Step 0 - transform the input level set into a source-type bump.
133  ParGridFunction source(&pfes);
134  source.ProjectCoefficient(zero_level_set);
135  // Optional smoothing of the initial level set.
137  // Transform so that the peak is at 0.
138  // Assumes range [-1, 1].
139  if (transform)
140  {
141  for (int i = 0; i < source.Size(); i++)
142  {
143  const double x = source(i);
144  source(i) = ((x < -1.0) || (x > 1.0)) ? 0.0 : (1.0 - x) * (1.0 + x);
145  }
146  }
147 
148  int amg_print_level = 0;
149 
150  // Solver.
151  CGSolver cg(MPI_COMM_WORLD);
152  cg.SetRelTol(1e-12);
153  cg.SetMaxIter(100);
155  OperatorPtr A;
156  Vector B, X;
157 
158  // Step 1 - diffuse.
159  ParGridFunction diffused_source(&pfes);
160  for (int i = 0; i < diffuse_iter; i++)
161  {
162  // Set up RHS.
163  ParLinearForm b(&pfes);
164  GridFunctionCoefficient src_coeff(&source);
165  b.AddDomainIntegrator(new DomainLFIntegrator(src_coeff));
166  b.Assemble();
167 
168  // Diffusion and mass terms in the LHS.
169  ParBilinearForm a_d(&pfes);
172  a_d.AddDomainIntegrator(new DiffusionIntegrator(t_coeff));
173  a_d.Assemble();
174 
175  // Solve with Dirichlet BC.
176  Array<int> ess_tdof_list;
177  if (pmesh.bdr_attributes.Size())
178  {
179  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
180  ess_bdr = 1;
181  pfes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
182  }
183  ParGridFunction u_dirichlet(&pfes);
184  u_dirichlet = 0.0;
185  a_d.FormLinearSystem(ess_tdof_list, u_dirichlet, b, A, X, B);
186  auto *prec = new HypreBoomerAMG;
187  prec->SetPrintLevel(amg_print_level);
188  cg.SetPreconditioner(*prec);
189  cg.SetOperator(*A);
190  cg.Mult(B, X);
191  a_d.RecoverFEMSolution(X, b, u_dirichlet);
192  delete prec;
193 
194  // Diffusion and mass terms in the LHS.
195  ParBilinearForm a_n(&pfes);
197  a_n.AddDomainIntegrator(new DiffusionIntegrator(t_coeff));
198  a_n.Assemble();
199 
200  // Solve with Neumann BC.
201  ParGridFunction u_neumann(&pfes);
202  ess_tdof_list.DeleteAll();
203  a_n.FormLinearSystem(ess_tdof_list, u_neumann, b, A, X, B);
204  auto *prec2 = new HypreBoomerAMG;
205  prec2->SetPrintLevel(amg_print_level);
206  cg.SetPreconditioner(*prec2);
207  cg.SetOperator(*A);
208  cg.Mult(B, X);
209  a_n.RecoverFEMSolution(X, b, u_neumann);
210  delete prec2;
211 
212  for (int ii = 0; ii < diffused_source.Size(); ii++)
213  {
214  // This assumes that the magnitudes of the two solutions are somewhat
215  // similar; otherwise one of the solutions would dominate and the BC
216  // won't look correct. To avoid this, it's good to have the source
217  // away from the boundary (i.e. have more resolution).
218  diffused_source(ii) = 0.5 * (u_neumann(ii) + u_dirichlet(ii));
219  }
220  source = diffused_source;
221  }
222 
223  // Step 2 - solve for the distance using the normalized gradient.
224  {
225  // RHS - normalized gradient.
226  ParLinearForm b2(&pfes);
227  NormalizedGradCoefficient grad_u(diffused_source, pmesh.Dimension());
229  b2.Assemble();
230 
231  // LHS - diffusion.
232  ParBilinearForm a2(&pfes);
234  a2.Assemble();
235 
236  // No BC.
237  Array<int> no_ess_tdofs;
238 
239  a2.FormLinearSystem(no_ess_tdofs, distance, b2, A, X, B);
240 
241  auto *prec = new HypreBoomerAMG;
242  prec->SetPrintLevel(amg_print_level);
243  OrthoSolver ortho(pfes.GetComm());
244  ortho.SetSolver(*prec);
245  cg.SetPreconditioner(ortho);
246  cg.SetOperator(*A);
247  cg.Mult(B, X);
248  a2.RecoverFEMSolution(X, b2, distance);
249  delete prec;
250  }
251 
252  // Shift the distance values to have minimum at zero.
253  double d_min_loc = distance.Min();
254  double d_min_glob;
255  MPI_Allreduce(&d_min_loc, &d_min_glob, 1, MPI_DOUBLE,
256  MPI_MIN, pfes.GetComm());
257  distance -= d_min_glob;
258 
259  if (vis_glvis)
260  {
261  char vishost[] = "localhost";
262  int visport = 19916;
263 
264  ParFiniteElementSpace fespace_vec(&pmesh, pfes.FEColl(),
265  pmesh.Dimension());
266  NormalizedGradCoefficient grad_u(diffused_source, pmesh.Dimension());
267  ParGridFunction x(&fespace_vec);
268  x.ProjectCoefficient(grad_u);
269 
270  socketstream sol_sock_x(vishost, visport);
271  sol_sock_x << "parallel " << pfes.GetNRanks() << " "
272  << pfes.GetMyRank() << "\n";
273  sol_sock_x.precision(8);
274  sol_sock_x << "solution\n" << pmesh << x;
275  sol_sock_x << "window_geometry " << 0 << " " << 0 << " "
276  << 500 << " " << 500 << "\n"
277  << "window_title '" << "Heat Directions" << "'\n"
278  << "keys evvRj*******A\n" << std::flush;
279  }
280 }
281 
282 double NormalizationDistanceSolver::NormalizationCoeff::
283 Eval(ElementTransformation &T,const IntegrationPoint &ip)
284 {
285  T.SetIntPoint(&ip);
286  Vector u_grad;
287  u.GetGradient(T, u_grad);
288  const double u_value = u.GetValue(T, ip);
289 
290  return u_value / sqrt(u_value * u_value + u_grad * u_grad + 1e-12);
291 }
292 
294  ParGridFunction &dist)
295 {
296  ParFiniteElementSpace &pfes = *dist.ParFESpace();
297 
298  ParGridFunction u_gf(&pfes);
299  u_gf.ProjectCoefficient(u_coeff);
300 
301  NormalizationCoeff rv_coeff(u_gf);
302  dist.ProjectDiscCoefficient(rv_coeff, GridFunction::AvgType::ARITHMETIC);
303 }
304 
306  ParGridFunction &fdist)
307 {
308  ParFiniteElementSpace* fesd=fdist.ParFESpace();
309 
310  auto check_h1 = dynamic_cast<const H1_FECollection *>(fesd->FEColl());
311  auto check_l2 = dynamic_cast<const L2_FECollection *>(fesd->FEColl());
312  MFEM_VERIFY((check_h1 || check_l2) && fesd->GetVDim() == 1,
313  "This solver supports only scalar H1 or L2 spaces.");
314 
315  ParMesh* mesh=fesd->GetParMesh();
316  const int dim=mesh->Dimension();
317 
318  MPI_Comm lcomm=fesd->GetComm();
319  int myrank;
320  MPI_Comm_rank(lcomm,&myrank);
321 
322  const int order = fesd->GetOrder(0);
323  H1_FECollection fecp(order, dim);
324  ParFiniteElementSpace fesp(mesh, &fecp, 1, Ordering::byVDIM);
325 
326  ParGridFunction wf(&fesp);
327  wf.ProjectCoefficient(func);
328  GradientGridFunctionCoefficient gf(&wf); //gradient of wf
329 
330 
331  ParGridFunction xf(&fesp);
332  HypreParVector *sv = xf.GetTrueDofs();
333  *sv=1.0;
334 
335  ParNonlinearForm *nf = new ParNonlinearForm(&fesp);
336 
337  PUMPLaplacian* pint = new PUMPLaplacian(&func,&gf,false);
338  nf->AddDomainIntegrator(pint);
339 
340  pint->SetPower(2);
341 
342  //define the solvers
343  HypreBoomerAMG *prec = new HypreBoomerAMG();
344  prec->SetPrintLevel(0);
345 
346  GMRESSolver *gmres;
347  gmres = new GMRESSolver(lcomm);
348  gmres->SetAbsTol(newton_abs_tol/10);
349  gmres->SetRelTol(newton_rel_tol/10);
350  gmres->SetMaxIter(100);
351  gmres->SetPrintLevel(0);
352  gmres->SetPreconditioner(*prec);
353 
354  NewtonSolver ns(lcomm);
355  ns.iterative_mode = true;
356  ns.SetSolver(*gmres);
357  ns.SetOperator(*nf);
359  ns.SetRelTol(newton_rel_tol);
360  ns.SetAbsTol(newton_abs_tol);
361  ns.SetMaxIter(newton_iter);
362 
363  Vector b; // RHS is zero
364  ns.Mult(b, *sv);
365 
366  for (int pp=3; pp<maxp; pp++)
367  {
368  if (myrank == 0 && (print_level.summary || print_level.iterations))
369  {
370  std::cout << "pp = " << pp << std::endl;
371  }
372  pint->SetPower(pp);
373  ns.Mult(b, *sv);
374  }
375 
376  xf.SetFromTrueDofs(*sv);
377  GridFunctionCoefficient gfx(&xf);
378  PProductCoefficient tsol(func,gfx);
379  fdist.ProjectCoefficient(tsol);
380 
381  delete gmres;
382  delete prec;
383  delete nf;
384  delete sv;
385 }
386 
387 
390  const Vector &elfun)
391 {
392  double energy = 0.0;
393  int ndof = el.GetDof();
394  int ndim = el.GetDim();
395  const IntegrationRule *ir = NULL;
396  int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
397  ir = &IntRules.Get(el.GetGeomType(), order);
398 
399  Vector shapef(ndof);
400  double fval;
401  double pval;
402  DenseMatrix B(ndof, ndim);
403  Vector qval(ndim);
404 
405  B=0.0;
406 
407  double w;
408  double ngrad2;
409 
410  for (int i = 0; i < ir->GetNPoints(); i++)
411  {
412  const IntegrationPoint &ip = ir->IntPoint(i);
413  trans.SetIntPoint(&ip);
414  w = trans.Weight();
415  w = ip.weight * w;
416 
417  fval=func->Eval(trans,ip);
418 
419  el.CalcPhysDShape(trans, B);
420  el.CalcPhysShape(trans,shapef);
421 
422  B.MultTranspose(elfun,qval);
423 
424  ngrad2=0.0;
425  for (int jj=0; jj<ndim; jj++)
426  {
427  ngrad2 = ngrad2 + qval(jj)*qval(jj);
428  }
429 
430  energy = energy + w * ngrad2 * diffcoef * 0.5;
431 
432  // add the external load -1 if fval > 0.0; 1 if fval < 0.0;
433  pval=shapef*elfun;
434 
435  energy = energy + w * pval * pval * 0.5;
436 
437  if (fval>0.0)
438  {
439  energy = energy - w*pval;
440  }
441  else if (fval<0.0)
442  {
443  energy = energy + w*pval;
444  }
445  }
446 
447  return energy;
448 }
449 
452  const Vector &elfun,
453  Vector &elvect)
454 {
455  int ndof = el.GetDof();
456  int ndim = el.GetDim();
457  const IntegrationRule *ir = NULL;
458  int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
459  ir = &IntRules.Get(el.GetGeomType(), order);
460 
461  elvect.SetSize(ndof);
462  elvect=0.0;
463 
464  Vector shapef(ndof);
465  double fval;
466  double pval;
467 
468  DenseMatrix B(ndof, ndim); //[diff_x,diff_y,diff_z]
469 
470  Vector qval(ndim); //[diff_x,diff_y,diff_z,u]
471  Vector lvec(ndof); //residual at ip
472 
473  B=0.0;
474  qval=0.0;
475 
476  double w;
477 
478  for (int i = 0; i < ir->GetNPoints(); i++)
479  {
480  const IntegrationPoint &ip = ir->IntPoint(i);
481  trans.SetIntPoint(&ip);
482  w = trans.Weight();
483  w = ip.weight * w;
484 
485  fval=func->Eval(trans,ip);
486 
487  el.CalcPhysDShape(trans, B);
488  el.CalcPhysShape(trans,shapef);
489 
490  B.MultTranspose(elfun,qval);
491  B.Mult(qval,lvec);
492 
493  elvect.Add(w * diffcoef,lvec);
494 
495  pval=shapef*elfun;
496 
497  elvect.Add(w * pval, shapef);
498 
499 
500  //add the load
501  //add the external load -1 if fval > 0.0; 1 if fval < 0.0;
502  pval=shapef*elfun;
503  if (fval>0.0)
504  {
505  elvect.Add( -w, shapef);
506  }
507  else if (fval<0.0)
508  {
509  elvect.Add( w, shapef);
510  }
511  }
512 }
513 
516  const Vector &elfun,
517  DenseMatrix &elmat)
518 {
519  int ndof = el.GetDof();
520  int ndim = el.GetDim();
521  const IntegrationRule *ir = NULL;
522  int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
523  ir = &IntRules.Get(el.GetGeomType(), order);
524 
525  elmat.SetSize(ndof,ndof);
526  elmat=0.0;
527 
528  Vector shapef(ndof);
529 
530  DenseMatrix B(ndof, ndim); //[diff_x,diff_y,diff_z]
531  B = 0.0;
532 
533  double w;
534 
535  for (int i = 0; i < ir->GetNPoints(); i++)
536  {
537  const IntegrationPoint &ip = ir->IntPoint(i);
538  trans.SetIntPoint(&ip);
539  w = trans.Weight();
540  w = ip.weight * w;
541 
542  el.CalcPhysDShape(trans, B);
543  el.CalcPhysShape(trans,shapef);
544 
545  AddMult_a_VVt(w, shapef, elmat);
546  AddMult_a_AAt(w * diffcoef, B, elmat);
547  }
548 }
549 
550 
553  const Vector &elfun)
554 {
555  double energy = 0.0;
556  int ndof = el.GetDof();
557  int ndim = el.GetDim();
558  const IntegrationRule *ir = NULL;
559  int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
560  ir = &IntRules.Get(el.GetGeomType(), order);
561 
562  Vector shapef(ndof);
563  double fval;
564  double pval;
565  double tval;
566  Vector vgrad(ndim);
567  DenseMatrix dshape(ndof, ndim);
568  DenseMatrix B(ndof, ndim);
569  Vector qval(ndim);
570  Vector tmpv(ndof);
571 
572  B=0.0;
573 
574  double w;
575  double ngrad2;
576 
577  for (int i = 0; i < ir->GetNPoints(); i++)
578  {
579  const IntegrationPoint &ip = ir->IntPoint(i);
580  trans.SetIntPoint(&ip);
581  w = trans.Weight();
582  w = ip.weight * w;
583 
584  fval=func->Eval(trans,ip);
585  fgrad->Eval(vgrad,trans,ip);
586  tval=fval;
587  if (fval<0.0)
588  {
589  fval=-fval;
590  vgrad*=-1.0;
591  }
592 
593  el.CalcPhysDShape(trans, dshape);
594  el.CalcPhysShape(trans,shapef);
595 
596  for (int jj=0; jj<ndim; jj++)
597  {
598  dshape.GetColumn(jj,tmpv);
599  tmpv*=fval;
600  tmpv.Add(vgrad[jj],shapef);
601  B.SetCol(jj,tmpv);
602  }
603  B.MultTranspose(elfun,qval);
604 
605  ngrad2=0.0;
606  for (int jj=0; jj<ndim; jj++)
607  {
608  ngrad2 = ngrad2 + qval(jj)*qval(jj);
609  }
610 
611  energy = energy + w * std::pow(ngrad2+ee*ee,pp/2.0)/pp;
612 
613  // add the external load -1 if fval > 0.0; 1 if fval < 0.0;
614  pval=shapef*elfun;
615  if (tval>0.0)
616  {
617  energy = energy - w * pval * tval;
618  }
619  else if (tval<0.0)
620  {
621  energy = energy + w * pval * tval;
622  }
623  }
624 
625  return energy;
626 }
627 
630  const Vector &elfun,
631  Vector &elvect)
632 {
633  int ndof = el.GetDof();
634  int ndim = el.GetDim();
635  const IntegrationRule *ir = NULL;
636  int order = 2 * el.GetOrder() + trans.OrderGrad(&el)+1;
637  ir = &IntRules.Get(el.GetGeomType(), order);
638 
639  elvect.SetSize(ndof);
640  elvect=0.0;
641 
642  Vector shapef(ndof);
643  double fval;
644  double tval;
645  Vector vgrad(3);
646 
647  DenseMatrix dshape(ndof, ndim);
648  DenseMatrix B(ndof, ndim); // [diff_x,diff_y,diff_z]
649 
650  Vector qval(ndim); // [diff_x,diff_y,diff_z,u]
651  Vector lvec(ndof); // residual at ip
652  Vector tmpv(ndof);
653 
654  B=0.0;
655  qval=0.0;
656 
657  double w;
658  double ngrad2;
659  double aa;
660 
661  for (int i = 0; i < ir->GetNPoints(); i++)
662  {
663  const IntegrationPoint &ip = ir->IntPoint(i);
664  trans.SetIntPoint(&ip);
665  w = trans.Weight();
666  w = ip.weight * w;
667 
668  fval=func->Eval(trans,ip);
669  fgrad->Eval(vgrad,trans,ip);
670  tval=fval;
671  if (fval<0.0)
672  {
673  fval=-fval;
674  vgrad*=-1.0;
675  }
676 
677  el.CalcPhysDShape(trans, dshape);
678  el.CalcPhysShape(trans,shapef);
679 
680  for (int jj=0; jj<ndim; jj++)
681  {
682  dshape.GetColumn(jj,tmpv);
683  tmpv*=fval;
684  tmpv.Add(vgrad[jj],shapef);
685  B.SetCol(jj,tmpv);
686  }
687 
688  B.MultTranspose(elfun,qval);
689 
690  ngrad2=0.0;
691  for (int jj=0; jj<ndim; jj++)
692  {
693  ngrad2 = ngrad2 + qval(jj)*qval(jj);
694  }
695 
696  aa = ngrad2 + ee*ee;
697  aa = std::pow(aa, (pp - 2.0) / 2.0);
698  B.Mult(qval,lvec);
699  elvect.Add(w * aa,lvec);
700 
701  // add the load
702  // add the external load -1 if tval > 0.0; 1 if tval < 0.0;
703  if (tval>0.0)
704  {
705  elvect.Add( -w*fval, shapef);
706  }
707  else if (tval<0.0)
708  {
709  elvect.Add( w*fval, shapef);
710  }
711  }
712 }
713 
716  const Vector &elfun,
717  DenseMatrix &elmat)
718 {
719  int ndof = el.GetDof();
720  int ndim = el.GetDim();
721  const IntegrationRule *ir = NULL;
722  int order = 2 * el.GetOrder() + trans.OrderGrad(&el)+1;
723  ir = &IntRules.Get(el.GetGeomType(), order);
724 
725  elmat.SetSize(ndof,ndof);
726  elmat=0.0;
727 
728  Vector shapef(ndof);
729  double fval;
730  Vector vgrad(ndim);
731 
732  Vector qval(ndim); // [diff_x,diff_y,diff_z,u]
733  DenseMatrix dshape(ndof, ndim);
734  DenseMatrix B(ndof, ndim); // [diff_x,diff_y,diff_z]
735  Vector lvec(ndof);
736  Vector tmpv(ndof);
737 
738  B=0.0;
739 
740  double w;
741  double ngrad2;
742  double aa, aa0, aa1;
743 
744  for (int i = 0; i < ir->GetNPoints(); i++)
745  {
746  const IntegrationPoint &ip = ir->IntPoint(i);
747  trans.SetIntPoint(&ip);
748  w = trans.Weight();
749  w = ip.weight * w;
750 
751  fval=func->Eval(trans,ip);
752  fgrad->Eval(vgrad,trans,ip);
753  if (fval<0.0)
754  {
755  fval=-fval;
756  vgrad*=-1.0;
757  }
758 
759  el.CalcPhysDShape(trans, dshape);
760  el.CalcPhysShape(trans,shapef);
761 
762  for (int jj=0; jj<ndim; jj++)
763  {
764  dshape.GetColumn(jj,tmpv);
765  tmpv*=fval;
766  tmpv.Add(vgrad[jj],shapef);
767  B.SetCol(jj,tmpv);
768  }
769 
770  B.MultTranspose(elfun,qval);
771  B.Mult(qval,lvec);
772 
773  ngrad2=0.0;
774  for (int jj=0; jj<ndim; jj++)
775  {
776  ngrad2 = ngrad2 + qval(jj)*qval(jj);
777  }
778 
779  aa = ngrad2 + ee * ee;
780  aa1 = std::pow(aa, (pp - 2.0) / 2.0);
781  aa0 = (pp-2.0) * std::pow(aa, (pp - 4.0) / 2.0);
782 
783  AddMult_a_VVt(w * aa0, lvec, elmat);
784  AddMult_a_AAt(w * aa1, B, elmat);
785  }
786 }
787 
788 
790 {
791  if (sint == nullptr)
792  {
793  sint = new ScreenedPoisson(func, rr);
794  nf->AddDomainIntegrator(sint);
795  *sv = 0.0;
796  gmres->SetOperator(nf->GetGradient(*sv));
797  }
798  else { sint->SetInput(func); }
799 
800  // form RHS
801  *sv = 0.0;
802  Vector rhs(sv->Size());
803  nf->Mult(*sv, rhs);
804  // filter the input field
805  gmres->Mult(rhs, *sv);
806 
807  gf.SetFromTrueDofs(*sv);
808  gf.Neg();
809 
810  GridFunctionCoefficient gfc(&gf);
811  ffield.ProjectCoefficient(gfc);
812 }
813 
814 }
Abstract class for all finite elements.
Definition: fe_base.hpp:232
const char vishost[]
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1032
Conjugate gradient method.
Definition: solvers.hpp:493
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void GetDerivative(int comp, int der_comp, ParGridFunction &der)
Parallel version of GridFunction::GetDerivative(); see its documentation.
Definition: pgridfunc.cpp:485
void SetCol(int c, const double *col)
Definition: densemat.cpp:2193
void ComputeScalarDistance(Coefficient &func, ParGridFunction &fdist)
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
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(...
Coefficient * func
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1108
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
int Dimension() const
Definition: mesh.hpp:1047
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
Class for domain integrator L(v) := (f, grad v)
Definition: lininteg.hpp:145
bool summary
A summary of the solver process will be reported after the last iteration to mfem::out.
Definition: solvers.hpp:91
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
virtual void Mult(const Vector &x, Vector &y) const
Evaluate the action of the NonlinearForm.
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:712
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:974
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Parallel non-linear operator on the true dofs.
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: solvers.hpp:689
virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
Project a discontinuous vector coefficient as a grid function on a continuous finite element space...
Definition: pgridfunc.cpp:546
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:525
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:686
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
Definition: hypre.cpp:3637
void DeleteAll()
Delete the whole array.
Definition: array.hpp:851
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
void ScalarDistToVector(ParGridFunction &dist_s, ParGridFunction &dist_v)
Definition: dist_solver.cpp:72
VectorCoefficient * fgrad
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:319
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect) override
Perform the local action of the NonlinearFormIntegrator.
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
void source(const Vector &x, Vector &f)
Definition: ex25.cpp:620
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
Class for parallel linear form.
Definition: plinearform.hpp:26
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:968
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:601
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:215
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
double b
Definition: lissajous.cpp:42
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1670
void SetInput(Coefficient &nfunc)
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
Vector coefficient defined as the Gradient of a scalar GridFunction.
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun) override
Compute the local energy.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:641
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
Definition: densemat.cpp:3233
Parallel smoothers in hypre.
Definition: hypre.hpp:971
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1331
const int visport
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:17
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
Definition: fe_base.cpp:191
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
IterativeSolver::PrintLevel print_level
Definition: dist_solver.hpp:29
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1817
void SetAbsTol(double atol)
Definition: solvers.hpp:200
void SetRelTol(double rtol)
Definition: solvers.hpp:199
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:310
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1194
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:571
virtual void ComputeVectorDistance(Coefficient &zero_level_set, ParGridFunction &distance)
GMRES method.
Definition: solvers.hpp:525
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:322
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect) override
Perform the local action of the NonlinearFormIntegrator.
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
void ComputeScalarDistance(Coefficient &zero_level_set, ParGridFunction &distance)
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
Definition: gridfunc.hpp:130
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
Solver wrapper which orthogonalizes the input and output vector.
Definition: solvers.hpp:1163
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
Definition: fe_base.cpp:181
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:253
Class for integration point with weight.
Definition: intrules.hpp:25
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:599
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Class for parallel bilinear form.
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun) override
Compute the local energy.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
Vector data type.
Definition: vector.hpp:60
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:143
void SetSolver(Solver &s)
Set the solver used by the OrthoSolver.
Definition: solvers.cpp:3431
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
bool iterations
Detailed information about each iteration will be reported to mfem::out.
Definition: solvers.hpp:88
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
void SetPower(double pp_)
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
void ComputeScalarDistance(Coefficient &u_coeff, ParGridFunction &dist)
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:326
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:3095
double GetElementVolume(int i)
Definition: mesh.cpp:112
void Filter(ParGridFunction &func, ParGridFunction &ffield)
virtual Operator & GetGradient(const Vector &x) const
Compute the gradient Operator of the NonlinearForm corresponding to the state x.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1806
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320
double AvgElementSize(ParMesh &pmesh)
Definition: dist_solver.cpp:42
void Neg()
(*this) = -(*this)
Definition: vector.cpp:306
virtual void ComputeScalarDistance(Coefficient &zero_level_set, ParGridFunction &distance)=0