1 // Copyright (c) 2010-2021, 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
DiffuseField(ParGridFunction & field,int smooth_steps)17 void DiffuseField(ParGridFunction &field, int smooth_steps)
18 {
19 // Setup the Laplacian operator.
20 ParBilinearForm *Lap = new ParBilinearForm(field.ParFESpace());
21 Lap->AddDomainIntegrator(new DiffusionIntegrator());
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
AvgElementSize(ParMesh & pmesh)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;
64 case Geometry::TETRAHEDRON:
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
ScalarDistToVector(ParGridFunction & dist_s,ParGridFunction & dist_v)72 void DistanceSolver::ScalarDistToVector(ParGridFunction &dist_s,
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
ComputeVectorDistance(Coefficient & zero_level_set,ParGridFunction & distance)107 void DistanceSolver::ComputeVectorDistance(Coefficient &zero_level_set,
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
ComputeScalarDistance(Coefficient & zero_level_set,ParGridFunction & distance)121 void HeatDistanceSolver::ComputeScalarDistance(Coefficient &zero_level_set,
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);
172 a_d.AddDomainIntegrator(new MassIntegrator);
173 ConstantCoefficient t_coeff(parameter_t);
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);
198 a_n.AddDomainIntegrator(new MassIntegrator);
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 i = 0; i < diffused_source.Size(); i++)
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(i) = 0.5 * (u_neumann(i) + u_dirichlet(i));
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());
230 b2.AddDomainIntegrator(new DomainLFGradIntegrator(grad_u));
231 b2.Assemble();
232
233 // LHS - diffusion.
234 ParBilinearForm a2(&pfes);
235 a2.AddDomainIntegrator(new DiffusionIntegrator);
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
ComputeScalarDistance(Coefficient & func,ParGridFunction & fdist)283 void PLapDistanceSolver::ComputeScalarDistance(Coefficient &func,
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
GetElementEnergy(const FiniteElement & el,ElementTransformation & trans,const Vector & elfun)369 double ScreenedPoisson::GetElementEnergy(const FiniteElement &el,
370 ElementTransformation &trans,
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
AssembleElementVector(const FiniteElement & el,ElementTransformation & trans,const Vector & elfun,Vector & elvect)431 void ScreenedPoisson::AssembleElementVector(const FiniteElement &el,
432 ElementTransformation &trans,
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
AssembleElementGrad(const FiniteElement & el,ElementTransformation & trans,const Vector & elfun,DenseMatrix & elmat)495 void ScreenedPoisson::AssembleElementGrad(const FiniteElement &el,
496 ElementTransformation &trans,
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
GetElementEnergy(const FiniteElement & el,ElementTransformation & trans,const Vector & elfun)532 double PUMPLaplacian::GetElementEnergy(const FiniteElement &el,
533 ElementTransformation &trans,
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
AssembleElementVector(const FiniteElement & el,ElementTransformation & trans,const Vector & elfun,Vector & elvect)609 void PUMPLaplacian::AssembleElementVector(const FiniteElement &el,
610 ElementTransformation &trans,
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
AssembleElementGrad(const FiniteElement & el,ElementTransformation & trans,const Vector & elfun,DenseMatrix & elmat)695 void PUMPLaplacian::AssembleElementGrad(const FiniteElement &el,
696 ElementTransformation &trans,
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
Filter(Coefficient & func,ParGridFunction & ffield)770 void PDEFilter::Filter(Coefficient &func, ParGridFunction &ffield)
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 }
796