MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
dist_solver.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 "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  const double eps = 1e-16;
97  for (int i = 0; i < size; i++)
98  {
99  const double vec_magn = std::sqrt(magn(i) + eps);
100  for (int d = 0; d < dim; d++)
101  {
102  dist_v(i + d*size) *= fabs(dist_s(i)) / vec_magn;
103  }
104  }
105 }
106 
108  ParGridFunction &distance)
109 {
110  ParFiniteElementSpace &pfes = *distance.ParFESpace();
111  MFEM_VERIFY(pfes.GetVDim() == pfes.GetMesh()->Dimension(),
112  "This function expects a vector ParGridFunction!");
113 
114  ParFiniteElementSpace pfes_s(pfes.GetParMesh(), pfes.FEColl());
115  ParGridFunction dist_s(&pfes_s);
116  ComputeScalarDistance(zero_level_set, dist_s);
117  ScalarDistToVector(dist_s, distance);
118 }
119 
120 
122  ParGridFunction &distance)
123 {
124  ParFiniteElementSpace &pfes = *distance.ParFESpace();
125 
126  auto check_h1 = dynamic_cast<const H1_FECollection *>(pfes.FEColl());
127  MFEM_VERIFY(check_h1 && pfes.GetVDim() == 1,
128  "This solver supports only scalar H1 spaces.");
129 
130  // Compute average mesh size (assumes similar cells).
131  ParMesh &pmesh = *pfes.GetParMesh();
132 
133  // Step 0 - transform the input level set into a source-type bump.
134  ParGridFunction source(&pfes);
135  source.ProjectCoefficient(zero_level_set);
136  // Optional smoothing of the initial level set.
137  if (smooth_steps > 0) { DiffuseField(source, smooth_steps); }
138  // Transform so that the peak is at 0.
139  // Assumes range [-1, 1].
140  if (transform)
141  {
142  for (int i = 0; i < source.Size(); i++)
143  {
144  const double x = source(i);
145  source(i) = ((x < -1.0) || (x > 1.0)) ? 0.0 : (1.0 - x) * (1.0 + x);
146  }
147  }
148 
149  const int cg_print_lvl = (print_level > 0) ? 1 : 0,
150  amg_print_lvl = (print_level > 1) ? 1 : 0;
151 
152  // Solver.
153  CGSolver cg(MPI_COMM_WORLD);
154  cg.SetRelTol(1e-12);
155  cg.SetMaxIter(100);
156  cg.SetPrintLevel(cg_print_lvl);
157  OperatorPtr A;
158  Vector B, X;
159 
160  // Step 1 - diffuse.
161  ParGridFunction diffused_source(&pfes);
162  for (int i = 0; i < diffuse_iter; i++)
163  {
164  // Set up RHS.
165  ParLinearForm b(&pfes);
166  GridFunctionCoefficient src_coeff(&source);
167  b.AddDomainIntegrator(new DomainLFIntegrator(src_coeff));
168  b.Assemble();
169 
170  // Diffusion and mass terms in the LHS.
171  ParBilinearForm a_d(&pfes);
174  a_d.AddDomainIntegrator(new DiffusionIntegrator(t_coeff));
175  a_d.Assemble();
176 
177  // Solve with Dirichlet BC.
178  Array<int> ess_tdof_list;
179  if (pmesh.bdr_attributes.Size())
180  {
181  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
182  ess_bdr = 1;
183  pfes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
184  }
185  ParGridFunction u_dirichlet(&pfes);
186  u_dirichlet = 0.0;
187  a_d.FormLinearSystem(ess_tdof_list, u_dirichlet, b, A, X, B);
188  auto *prec = new HypreBoomerAMG;
189  prec->SetPrintLevel(amg_print_lvl);
190  cg.SetPreconditioner(*prec);
191  cg.SetOperator(*A);
192  cg.Mult(B, X);
193  a_d.RecoverFEMSolution(X, b, u_dirichlet);
194  delete prec;
195 
196  // Diffusion and mass terms in the LHS.
197  ParBilinearForm a_n(&pfes);
199  a_n.AddDomainIntegrator(new DiffusionIntegrator(t_coeff));
200  a_n.Assemble();
201 
202  // Solve with Neumann BC.
203  ParGridFunction u_neumann(&pfes);
204  ess_tdof_list.DeleteAll();
205  a_n.FormLinearSystem(ess_tdof_list, u_neumann, b, A, X, B);
206  auto *prec2 = new HypreBoomerAMG;
207  prec2->SetPrintLevel(amg_print_lvl);
208  cg.SetPreconditioner(*prec2);
209  cg.SetOperator(*A);
210  cg.Mult(B, X);
211  a_n.RecoverFEMSolution(X, b, u_neumann);
212  delete prec2;
213 
214  for (int ii = 0; ii < diffused_source.Size(); ii++)
215  {
216  // This assumes that the magnitudes of the two solutions are somewhat
217  // similar; otherwise one of the solutions would dominate and the BC
218  // won't look correct. To avoid this, it's good to have the source
219  // away from the boundary (i.e. have more resolution).
220  diffused_source(ii) = 0.5 * (u_neumann(ii) + u_dirichlet(ii));
221  }
222  source = diffused_source;
223  }
224 
225  // Step 2 - solve for the distance using the normalized gradient.
226  {
227  // RHS - normalized gradient.
228  ParLinearForm b2(&pfes);
229  NormalizedGradCoefficient grad_u(diffused_source, pmesh.Dimension());
231  b2.Assemble();
232 
233  // LHS - diffusion.
234  ParBilinearForm a2(&pfes);
236  a2.Assemble();
237 
238  // No BC.
239  Array<int> no_ess_tdofs;
240 
241  a2.FormLinearSystem(no_ess_tdofs, distance, b2, A, X, B);
242 
243  auto *prec = new HypreBoomerAMG;
244  prec->SetPrintLevel(amg_print_lvl);
245  cg.SetPreconditioner(*prec);
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 
284  ParGridFunction &fdist)
285 {
286  ParFiniteElementSpace* fesd=fdist.ParFESpace();
287 
288  auto check_h1 = dynamic_cast<const H1_FECollection *>(fesd->FEColl());
289  auto check_l2 = dynamic_cast<const L2_FECollection *>(fesd->FEColl());
290  MFEM_VERIFY((check_h1 || check_l2) && fesd->GetVDim() == 1,
291  "This solver supports only scalar H1 or L2 spaces.");
292 
293  ParMesh* mesh=fesd->GetParMesh();
294  const int dim=mesh->Dimension();
295 
296  MPI_Comm lcomm=fesd->GetComm();
297  int myrank;
298  MPI_Comm_rank(lcomm,&myrank);
299 
300  const int order = fesd->GetOrder(0);
301  H1_FECollection fecp(order, dim);
302  ParFiniteElementSpace fesp(mesh, &fecp, 1, Ordering::byVDIM);
303 
304  ParGridFunction wf(&fesp);
305  wf.ProjectCoefficient(func);
306  GradientGridFunctionCoefficient gf(&wf); //gradient of wf
307 
308 
309  ParGridFunction xf(&fesp);
310  HypreParVector *sv = xf.GetTrueDofs();
311  *sv=1.0;
312 
313  ParNonlinearForm *nf = new ParNonlinearForm(&fesp);
314 
315  PUMPLaplacian* pint = new PUMPLaplacian(&func,&gf,false);
316  nf->AddDomainIntegrator(pint);
317 
318  pint->SetPower(2);
319 
320  //define the solvers
321  HypreBoomerAMG *prec = new HypreBoomerAMG();
322  prec->SetPrintLevel((print_level > 1) ? 1 : 0);
323 
324  GMRESSolver *gmres;
325  gmres = new GMRESSolver(lcomm);
326  gmres->SetAbsTol(newton_abs_tol/10);
327  gmres->SetRelTol(newton_rel_tol/10);
328  gmres->SetMaxIter(100);
329  gmres->SetPrintLevel((print_level > 1) ? 1 : 0);
330  gmres->SetPreconditioner(*prec);
331 
332  NewtonSolver ns(lcomm);
333  ns.iterative_mode = true;
334  ns.SetSolver(*gmres);
335  ns.SetOperator(*nf);
336  ns.SetPrintLevel((print_level == 0) ? -1 : 0);
337  ns.SetRelTol(newton_rel_tol);
338  ns.SetAbsTol(newton_abs_tol);
339  ns.SetMaxIter(newton_iter);
340 
341  Vector b; // RHS is zero
342  ns.Mult(b, *sv);
343 
344  for (int pp=3; pp<maxp; pp++)
345  {
346  if (myrank == 0 && print_level > 0) { std::cout<<"pp="<<pp<<std::endl; }
347  pint->SetPower(pp);
348  ns.Mult(b, *sv);
349  }
350 
351  xf.SetFromTrueDofs(*sv);
352  GridFunctionCoefficient gfx(&xf);
353  PProductCoefficient tsol(func,gfx);
354  fdist.ProjectCoefficient(tsol);
355 
356  // (optional) Force positive distances everywhere.
357  // for (int i = 0; i < fdist.Size(); i++)
358  // {
359  // fdist(i) = fabs(fdist(i));
360  // }
361 
362  delete gmres;
363  delete prec;
364  delete nf;
365  delete sv;
366 }
367 
368 
371  const Vector &elfun)
372 {
373  double energy = 0.0;
374  int ndof = el.GetDof();
375  int ndim = el.GetDim();
376  const IntegrationRule *ir = NULL;
377  int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
378  ir = &IntRules.Get(el.GetGeomType(), order);
379 
380  Vector shapef(ndof);
381  double fval;
382  double pval;
383  DenseMatrix B(ndof, ndim);
384  Vector qval(ndim);
385 
386  B=0.0;
387 
388  double w;
389  double ngrad2;
390 
391  for (int i = 0; i < ir->GetNPoints(); i++)
392  {
393  const IntegrationPoint &ip = ir->IntPoint(i);
394  trans.SetIntPoint(&ip);
395  w = trans.Weight();
396  w = ip.weight * w;
397 
398  fval=func->Eval(trans,ip);
399 
400  el.CalcPhysDShape(trans, B);
401  el.CalcPhysShape(trans,shapef);
402 
403  B.MultTranspose(elfun,qval);
404 
405  ngrad2=0.0;
406  for (int jj=0; jj<ndim; jj++)
407  {
408  ngrad2 = ngrad2 + qval(jj)*qval(jj);
409  }
410 
411  energy = energy + w * ngrad2 * diffcoef * 0.5;
412 
413  // add the external load -1 if fval > 0.0; 1 if fval < 0.0;
414  pval=shapef*elfun;
415 
416  energy = energy + w * pval * pval * 0.5;
417 
418  if (fval>0.0)
419  {
420  energy = energy - w*pval;
421  }
422  else if (fval<0.0)
423  {
424  energy = energy + w*pval;
425  }
426  }
427 
428  return energy;
429 }
430 
433  const Vector &elfun,
434  Vector &elvect)
435 {
436  int ndof = el.GetDof();
437  int ndim = el.GetDim();
438  const IntegrationRule *ir = NULL;
439  int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
440  ir = &IntRules.Get(el.GetGeomType(), order);
441 
442  elvect.SetSize(ndof);
443  elvect=0.0;
444 
445  Vector shapef(ndof);
446  double fval;
447  double pval;
448 
449  DenseMatrix B(ndof, ndim); //[diff_x,diff_y,diff_z]
450 
451  Vector qval(ndim); //[diff_x,diff_y,diff_z,u]
452  Vector lvec(ndof); //residual at ip
453 
454  B=0.0;
455  qval=0.0;
456 
457  double w;
458 
459  for (int i = 0; i < ir->GetNPoints(); i++)
460  {
461  const IntegrationPoint &ip = ir->IntPoint(i);
462  trans.SetIntPoint(&ip);
463  w = trans.Weight();
464  w = ip.weight * w;
465 
466  fval=func->Eval(trans,ip);
467 
468  el.CalcPhysDShape(trans, B);
469  el.CalcPhysShape(trans,shapef);
470 
471  B.MultTranspose(elfun,qval);
472  B.Mult(qval,lvec);
473 
474  elvect.Add(w * diffcoef,lvec);
475 
476  pval=shapef*elfun;
477 
478  elvect.Add(w * pval, shapef);
479 
480 
481  //add the load
482  //add the external load -1 if fval > 0.0; 1 if fval < 0.0;
483  pval=shapef*elfun;
484  if (fval>0.0)
485  {
486  elvect.Add( -w, shapef);
487  }
488  else if (fval<0.0)
489  {
490  elvect.Add( w, shapef);
491  }
492  }
493 }
494 
497  const Vector &elfun,
498  DenseMatrix &elmat)
499 {
500  int ndof = el.GetDof();
501  int ndim = el.GetDim();
502  const IntegrationRule *ir = NULL;
503  int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
504  ir = &IntRules.Get(el.GetGeomType(), order);
505 
506  elmat.SetSize(ndof,ndof);
507  elmat=0.0;
508 
509  Vector shapef(ndof);
510 
511  DenseMatrix B(ndof, ndim); //[diff_x,diff_y,diff_z]
512  B = 0.0;
513 
514  double w;
515 
516  for (int i = 0; i < ir->GetNPoints(); i++)
517  {
518  const IntegrationPoint &ip = ir->IntPoint(i);
519  trans.SetIntPoint(&ip);
520  w = trans.Weight();
521  w = ip.weight * w;
522 
523  el.CalcPhysDShape(trans, B);
524  el.CalcPhysShape(trans,shapef);
525 
526  AddMult_a_VVt(w, shapef, elmat);
527  AddMult_a_AAt(w * diffcoef, B, elmat);
528  }
529 }
530 
531 
534  const Vector &elfun)
535 {
536  double energy = 0.0;
537  int ndof = el.GetDof();
538  int ndim = el.GetDim();
539  const IntegrationRule *ir = NULL;
540  int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
541  ir = &IntRules.Get(el.GetGeomType(), order);
542 
543  Vector shapef(ndof);
544  double fval;
545  double pval;
546  double tval;
547  Vector vgrad(ndim);
548  DenseMatrix dshape(ndof, ndim);
549  DenseMatrix B(ndof, ndim);
550  Vector qval(ndim);
551  Vector tmpv(ndof);
552 
553  B=0.0;
554 
555  double w;
556  double ngrad2;
557 
558  for (int i = 0; i < ir->GetNPoints(); i++)
559  {
560  const IntegrationPoint &ip = ir->IntPoint(i);
561  trans.SetIntPoint(&ip);
562  w = trans.Weight();
563  w = ip.weight * w;
564 
565  fval=func->Eval(trans,ip);
566  fgrad->Eval(vgrad,trans,ip);
567  tval=fval;
568  if (fval<0.0)
569  {
570  fval=-fval;
571  vgrad*=-1.0;
572  }
573 
574  el.CalcPhysDShape(trans, dshape);
575  el.CalcPhysShape(trans,shapef);
576 
577  for (int jj=0; jj<ndim; jj++)
578  {
579  dshape.GetColumn(jj,tmpv);
580  tmpv*=fval;
581  tmpv.Add(vgrad[jj],shapef);
582  B.SetCol(jj,tmpv);
583  }
584  B.MultTranspose(elfun,qval);
585 
586  ngrad2=0.0;
587  for (int jj=0; jj<ndim; jj++)
588  {
589  ngrad2 = ngrad2 + qval(jj)*qval(jj);
590  }
591 
592  energy = energy + w * std::pow(ngrad2+ee*ee,pp/2.0)/pp;
593 
594  // add the external load -1 if fval > 0.0; 1 if fval < 0.0;
595  pval=shapef*elfun;
596  if (tval>0.0)
597  {
598  energy = energy - w * pval * tval;
599  }
600  else if (tval<0.0)
601  {
602  energy = energy + w * pval * tval;
603  }
604  }
605 
606  return energy;
607 }
608 
611  const Vector &elfun,
612  Vector &elvect)
613 {
614  int ndof = el.GetDof();
615  int ndim = el.GetDim();
616  const IntegrationRule *ir = NULL;
617  int order = 2 * el.GetOrder() + trans.OrderGrad(&el)+1;
618  ir = &IntRules.Get(el.GetGeomType(), order);
619 
620  elvect.SetSize(ndof);
621  elvect=0.0;
622 
623  Vector shapef(ndof);
624  double fval;
625  double tval;
626  Vector vgrad(3);
627 
628  DenseMatrix dshape(ndof, ndim);
629  DenseMatrix B(ndof, ndim); // [diff_x,diff_y,diff_z]
630 
631  Vector qval(ndim); // [diff_x,diff_y,diff_z,u]
632  Vector lvec(ndof); // residual at ip
633  Vector tmpv(ndof);
634 
635  B=0.0;
636  qval=0.0;
637 
638  double w;
639  double ngrad2;
640  double aa;
641 
642  for (int i = 0; i < ir->GetNPoints(); i++)
643  {
644  const IntegrationPoint &ip = ir->IntPoint(i);
645  trans.SetIntPoint(&ip);
646  w = trans.Weight();
647  w = ip.weight * w;
648 
649  fval=func->Eval(trans,ip);
650  fgrad->Eval(vgrad,trans,ip);
651  tval=fval;
652  if (fval<0.0)
653  {
654  fval=-fval;
655  vgrad*=-1.0;
656  }
657 
658  el.CalcPhysDShape(trans, dshape);
659  el.CalcPhysShape(trans,shapef);
660 
661  for (int jj=0; jj<ndim; jj++)
662  {
663  dshape.GetColumn(jj,tmpv);
664  tmpv*=fval;
665  tmpv.Add(vgrad[jj],shapef);
666  B.SetCol(jj,tmpv);
667  }
668 
669  B.MultTranspose(elfun,qval);
670 
671  ngrad2=0.0;
672  for (int jj=0; jj<ndim; jj++)
673  {
674  ngrad2 = ngrad2 + qval(jj)*qval(jj);
675  }
676 
677  aa = ngrad2 + ee*ee;
678  aa = std::pow(aa, (pp - 2.0) / 2.0);
679  B.Mult(qval,lvec);
680  elvect.Add(w * aa,lvec);
681 
682  // add the load
683  // add the external load -1 if tval > 0.0; 1 if tval < 0.0;
684  if (tval>0.0)
685  {
686  elvect.Add( -w*fval, shapef);
687  }
688  else if (tval<0.0)
689  {
690  elvect.Add( w*fval, shapef);
691  }
692  }
693 }
694 
697  const Vector &elfun,
698  DenseMatrix &elmat)
699 {
700  int ndof = el.GetDof();
701  int ndim = el.GetDim();
702  const IntegrationRule *ir = NULL;
703  int order = 2 * el.GetOrder() + trans.OrderGrad(&el)+1;
704  ir = &IntRules.Get(el.GetGeomType(), order);
705 
706  elmat.SetSize(ndof,ndof);
707  elmat=0.0;
708 
709  Vector shapef(ndof);
710  double fval;
711  Vector vgrad(ndim);
712 
713  Vector qval(ndim); // [diff_x,diff_y,diff_z,u]
714  DenseMatrix dshape(ndof, ndim);
715  DenseMatrix B(ndof, ndim); // [diff_x,diff_y,diff_z]
716  Vector lvec(ndof);
717  Vector tmpv(ndof);
718 
719  B=0.0;
720 
721  double w;
722  double ngrad2;
723  double aa, aa0, aa1;
724 
725  for (int i = 0; i < ir->GetNPoints(); i++)
726  {
727  const IntegrationPoint &ip = ir->IntPoint(i);
728  trans.SetIntPoint(&ip);
729  w = trans.Weight();
730  w = ip.weight * w;
731 
732  fval=func->Eval(trans,ip);
733  fgrad->Eval(vgrad,trans,ip);
734  if (fval<0.0)
735  {
736  fval=-fval;
737  vgrad*=-1.0;
738  }
739 
740  el.CalcPhysDShape(trans, dshape);
741  el.CalcPhysShape(trans,shapef);
742 
743  for (int jj=0; jj<ndim; jj++)
744  {
745  dshape.GetColumn(jj,tmpv);
746  tmpv*=fval;
747  tmpv.Add(vgrad[jj],shapef);
748  B.SetCol(jj,tmpv);
749  }
750 
751  B.MultTranspose(elfun,qval);
752  B.Mult(qval,lvec);
753 
754  ngrad2=0.0;
755  for (int jj=0; jj<ndim; jj++)
756  {
757  ngrad2 = ngrad2 + qval(jj)*qval(jj);
758  }
759 
760  aa = ngrad2 + ee * ee;
761  aa1 = std::pow(aa, (pp - 2.0) / 2.0);
762  aa0 = (pp-2.0) * std::pow(aa, (pp - 4.0) / 2.0);
763 
764  AddMult_a_VVt(w * aa0, lvec, elmat);
765  AddMult_a_AAt(w * aa1, B, elmat);
766  }
767 }
768 
769 
771 {
772  if (sint == nullptr)
773  {
774  sint = new ScreenedPoisson(func, rr);
775  nf->AddDomainIntegrator(sint);
776  *sv = 0.0;
777  gmres->SetOperator(nf->GetGradient(*sv));
778  }
779  else { sint->SetInput(func); }
780 
781  // form RHS
782  *sv = 0.0;
783  Vector rhs(sv->Size());
784  nf->Mult(*sv, rhs);
785  // filter the input field
786  gmres->Mult(rhs, *sv);
787 
788  gf.SetFromTrueDofs(*sv);
789  gf.Neg();
790 
791  GridFunctionCoefficient gfc(&gf);
792  ffield.ProjectCoefficient(gfc);
793 }
794 
795 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe_base.hpp:235
const char vishost[]
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:575
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
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:465
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual int OrderGrad(const FiniteElement *fe) const =0
Return the order of .
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:1788
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
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:960
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
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:711
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 ...
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
Class for domain integrator L(v) := (f, grad v)
Definition: lininteg.hpp:127
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Parallel non-linear operator on the true dofs.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: solvers.hpp:661
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:655
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1062
void DeleteAll()
Delete the whole array.
Definition: array.hpp:846
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
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:320
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect) override
Perform the local action of the NonlinearFormIntegrator.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
Definition: hypre.cpp:3426
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:971
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1443
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 MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:188
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:433
double b
Definition: lissajous.cpp:42
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1523
void SetInput(Coefficient &nfunc)
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
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.
virtual void Mult(const Vector &x, Vector &y) const
Evaluate the action of the NonlinearForm.
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:613
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
Definition: densemat.cpp:2828
Parallel smoothers in hypre.
Definition: hypre.hpp:896
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
const int visport
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:557
int Dimension() const
Definition: mesh.hpp:999
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:17
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1270
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1804
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
Definition: gridfunc.hpp:134
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:148
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1182
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
void SetAbsTol(double atol)
Definition: solvers.hpp:199
void SetRelTol(double rtol)
Definition: solvers.hpp:198
MPI_Comm GetComm() const
Definition: pmesh.hpp:288
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
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
virtual void ComputeVectorDistance(Coefficient &zero_level_set, ParGridFunction &distance)
int GetOrder(int i) const
Returns the polynomial degree of the i&#39;th finite element.
Definition: fespace.hpp:547
GMRES method.
Definition: solvers.hpp:497
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)
virtual Operator & GetGradient(const Vector &x) const
Compute the gradient Operator of the NonlinearForm corresponding to the state x.
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:252
Class for integration point with weight.
Definition: intrules.hpp:25
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:143
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.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
Vector data type.
Definition: vector.hpp:60
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:577
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:160
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:216
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:337
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
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:2690
double GetElementVolume(int i)
Definition: mesh.cpp:112
void Filter(ParGridFunction &func, ParGridFunction &ffield)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1793
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:284
double AvgElementSize(ParMesh &pmesh)
Definition: dist_solver.cpp:42
void Neg()
(*this) = -(*this)
Definition: vector.cpp:292
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
virtual void ComputeScalarDistance(Coefficient &zero_level_set, ParGridFunction &distance)=0