MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
sbm_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 "sbm_solver.hpp"
13 
14 namespace mfem
15 {
16 
18  const IntegrationPoint & ip,
19  const Vector &D)
20 {
21  if (constantcoefficient) { return constant; }
22 
23  Vector transip;
24  T.Transform(ip, transip);
25  for (int i = 0; i < D.Size(); i++)
26  {
27  transip(i) += D(i);
28  }
29 
30  return Function(transip);
31 }
32 
35  const IntegrationPoint & ip,
36  const Vector &D)
37 {
38  Vector transip;
39  T.Transform(ip, transip);
40  for (int i = 0; i < D.Size(); i++)
41  {
42  transip(i) += D(i);
43  }
44 
45  Function(transip, V);
46 }
47 
49  const FiniteElement &el1, const FiniteElement &el2,
51 {
52  int dim, ndof1, ndof2, ndof, ndoftotal;
53  double w;
54  DenseMatrix temp_elmat;
55 
56  dim = el1.GetDim();
57  ndof1 = el1.GetDof();
58  ndof2 = el2.GetDof();
59  ndoftotal = Trans.ElementType == ElementTransformation::BDR_FACE ?
60  ndof1 : ndof1 + ndof2;
61 
62  elmat.SetSize(ndoftotal);
63  elmat = 0.0;
64 
65  bool elem1f = true; // flag indicating whether Trans.Elem1No is part of the
66  // surrogate domain or not.
67  int elem1 = Trans.Elem1No,
68  elem2 = Trans.Elem2No,
69  marker1 = (*elem_marker)[elem1];
70 
71  int marker2;
72  if (Trans.Elem2No >= NEproc)
73  {
74  marker2 = (*elem_marker)[NEproc+par_shared_face_count];
76  }
78  {
79  marker2 = marker1;
80  }
81  else
82  {
83  marker2 = (*elem_marker)[elem2];
84  }
85 
86  if (!include_cut_cell)
87  {
88  // 1 is inside and 2 is cut or 1 is a boundary element.
89  if (marker1 == ShiftedFaceMarker::SBElementType::INSIDE &&
90  (cut_marker.Find(marker2) != -1 ||
92  {
93  elem1f = true;
94  }
95  // 1 is cut, 2 is inside
96  else if (cut_marker.Find(marker1) != -1 &&
97  marker2 == ShiftedFaceMarker::SBElementType::INSIDE)
98  {
99  if (Trans.Elem2No >= NEproc) { return; }
100  elem1f = false;
101  }
102  else
103  {
104  return;
105  }
106  }
107  else
108  {
109  // 1 is cut and 2 is outside or 1 is a boundary element.
110  if (cut_marker.Find(marker1) != -1 &&
111  (marker2 == ShiftedFaceMarker::SBElementType::OUTSIDE ||
113  {
114  elem1f = true;
115  }
116  // 1 is outside, 2 is cut
117  else if (marker1 == ShiftedFaceMarker::SBElementType::OUTSIDE &&
118  cut_marker.Find(marker2) != -1)
119  {
120  if (Trans.Elem2No >= NEproc) { return; }
121  elem1f = false;
122  }
123  else
124  {
125  return;
126  }
127  }
128 
129  ndof = elem1f ? ndof1 : ndof2;
130 
131  temp_elmat.SetSize(ndof);
132  temp_elmat = 0.;
133 
134  nor.SetSize(dim);
135  nh.SetSize(dim);
136  ni.SetSize(dim);
137  adjJ.SetSize(dim);
138 
139  shape.SetSize(ndof);
140  dshape.SetSize(ndof, dim);
141  dshapephys.SetSize(ndof, dim);
142  Vector dshapephysdd(ndof);
143  dshapedn.SetSize(ndof);
144  Vector wrk = shape;
145 
146  const IntegrationRule *ir = IntRule;
147  if (ir == NULL)
148  {
149  int order = elem1f ? 4*el1.GetOrder() : 4*el2.GetOrder();
150  ir = &IntRules.Get(Trans.GetGeometryType(), order);
151  }
152 
153  Array<DenseMatrix *> dkphi_dxk;
154  DenseMatrix grad_phys;
155  Vector Factorial;
156  Array<DenseMatrix *> grad_phys_dir;
157 
158  if (nterms > 0)
159  {
160  if (elem1f)
161  {
162  el1.ProjectGrad(el1, *Trans.Elem1, grad_phys);
163  }
164  else
165  {
166  el2.ProjectGrad(el2, *Trans.Elem2, grad_phys);
167  }
168 
169  DenseMatrix grad_work;
170  grad_phys_dir.SetSize(dim); // NxN matrices for derivative in each direction
171  for (int i = 0; i < dim; i++)
172  {
173  grad_phys_dir[i] = new DenseMatrix(ndof, ndof);
174  grad_phys_dir[i]->CopyRows(grad_phys, i*ndof, (i+1)*ndof-1);
175  }
176 
177  DenseMatrix grad_phys_work = grad_phys;
178  grad_phys_work.SetSize(ndof, ndof*dim);
179 
180  dkphi_dxk.SetSize(nterms);
181 
182  for (int i = 0; i < nterms; i++)
183  {
184  int sz1 = pow(dim, i+1);
185  dkphi_dxk[i] = new DenseMatrix(ndof, ndof*sz1*dim);
186  int loc_col_per_dof = sz1;
187  int tot_col_per_dof = loc_col_per_dof*dim;
188  for (int k = 0; k < dim; k++)
189  {
190  grad_work.SetSize(ndof, ndof*sz1);
191  // grad_work[k] has derivative in kth direction for each DOF.
192  // grad_work[0] has d^2phi/dx^2 and d^2phi/dxdy terms and
193  // grad_work[1] has d^2phi/dydx and d^2phi/dy2 terms for each dof
194  if (i == 0)
195  {
196  Mult(*grad_phys_dir[k], grad_phys_work, grad_work);
197  }
198  else
199  {
200  Mult(*grad_phys_dir[k], *dkphi_dxk[i-1], grad_work);
201  }
202  // Now we must place columns for each dof together so that they are
203  // in order: d^2phi/dx^2, d^2phi/dxdy, d^2phi/dydx, d^2phi/dy2.
204  for (int j = 0; j < ndof; j++)
205  {
206  for (int d = 0; d < loc_col_per_dof; d++)
207  {
208  Vector col;
209  grad_work.GetColumn(j*loc_col_per_dof+d, col);
210  dkphi_dxk[i]->SetCol(j*tot_col_per_dof+k*loc_col_per_dof+d, col);
211  }
212  }
213  }
214  }
215 
216  for (int i = 0; i < grad_phys_dir.Size(); i++)
217  {
218  delete grad_phys_dir[i];
219  }
220 
221  Factorial.SetSize(nterms);
222  Factorial(0) = 2;
223  for (int i = 1; i < nterms; i++)
224  {
225  Factorial(i) = Factorial(i-1)*(i+2);
226  }
227  }
228 
229  DenseMatrix q_hess_dn(dim, ndof);
230  Vector q_hess_dn_work(q_hess_dn.GetData(), ndof*dim);
231  Vector q_hess_dot_d(ndof);
232 
233  Vector D(vD->GetVDim());
234  // Assemble: -< \nabla u.n, w >
235  // -< u + \nabla u.d + h.o.t, \nabla w.n>
236  // +<alpha h^{-1} (u + \nabla u.d + h.o.t), w + \nabla w.d + h.o.t>
237  for (int p = 0; p < ir->GetNPoints(); p++)
238  {
239  const IntegrationPoint &ip = ir->IntPoint(p);
240 
241  // Set the integration point in the face and the neighboring elements
242  Trans.SetAllIntPoints(&ip);
243 
244  // Access the neighboring elements' integration points
245  // Note: eip2 will only contain valid data if Elem2 exists
246  const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
247  const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
248 
249  if (dim == 1)
250  {
251  nor(0) = 2*eip1.x - 1.0;
252  }
253  else
254  {
255  // Note: this normal accounts for the weight of the surface transformation
256  // Jacobian i.e. nor = nhat*det(J)
257  CalcOrtho(Trans.Jacobian(), nor);
258  }
259  vD->Eval(D, Trans, ip);
260 
261  // Make sure the normal vector is pointing outside the domain.
262  if (!elem1f) { nor *= -1; }
263 
264  if (elem1f)
265  {
266  el1.CalcShape(eip1, shape);
267  el1.CalcDShape(eip1, dshape);
268  w = ip.weight/Trans.Elem1->Weight();
269  CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
270  }
271  else
272  {
273  el1.CalcShape(eip2, shape);
274  el1.CalcDShape(eip2, dshape);
275  w = ip.weight/Trans.Elem2->Weight();
276  CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
277  }
278 
279  ni.Set(w, nor); // alpha_k*nor/det(J)
280  adjJ.Mult(ni, nh);
281  dshape.Mult(nh, dshapedn); // dphi/dn * Jinv * alpha_k * nor
282 
283  // <grad u.n, w> - Term 2
284  AddMult_a_VWt(-1., shape, dshapedn, temp_elmat);
285 
286  if (elem1f) { el1.CalcPhysDShape(*(Trans.Elem1), dshapephys); }
287  else { el1.CalcPhysDShape(*(Trans.Elem2), dshapephys); } // dphi/dx
288  dshapephys.Mult(D, dshapephysdd); // dphi/dx.D);
289 
290  q_hess_dot_d = 0.;
291  for (int i = 0; i < nterms; i++)
292  {
293  int sz1 = pow(dim, i+1);
294  DenseMatrix T1(dim, ndof*sz1);
295  Vector T1_wrk(T1.GetData(), dim*ndof*sz1);
296  dkphi_dxk[i]->MultTranspose(shape, T1_wrk);
297 
298  DenseMatrix T2;
299  Vector T2_wrk;
300  for (int j = 0; j < i+1; j++)
301  {
302  int sz2 = pow(dim, i-j);
303  T2.SetSize(dim, ndof*sz2);
304  T2_wrk.SetDataAndSize(T2.GetData(), dim*ndof*sz2);
305  T1.MultTranspose(D, T2_wrk);
306  T1 = T2;
307  }
308  Vector q_hess_dot_d_work(ndof);
309  T1.MultTranspose(D, q_hess_dot_d_work);
310  q_hess_dot_d_work *= 1./Factorial(i);
311  q_hess_dot_d += q_hess_dot_d_work;
312  }
313 
314  wrk = shape;
315  wrk += dshapephysdd;
316  wrk += q_hess_dot_d;
317  // <u + grad u.d + h.o.t, grad w.n> - Term 3
318  AddMult_a_VWt(-1., dshapedn, wrk, temp_elmat);
319 
320  double hinvdx;
321  if (elem1f) { hinvdx = nor*nor/Trans.Elem1->Weight(); }
322  else { hinvdx = nor*nor/Trans.Elem2->Weight(); }
323 
324  w = ip.weight*alpha*hinvdx;
325  // + <alpha * hinv * u + grad u.d + h.o.t, w + grad w.d + h.o.t> - Term 4
326  AddMult_a_VVt(w, wrk, temp_elmat);
327  } // p < ir->GetNPoints()
328  int offset = elem1f ? 0 : ndof1;
329  elmat.CopyMN(temp_elmat, offset, offset);
330 
331  for (int i = 0; i < dkphi_dxk.Size(); i++)
332  {
333  delete dkphi_dxk[i];
334  }
335 }
336 
337 
339  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
340 {
341  mfem_error("SBM2DirichletLFIntegrator::AssembleRHSElementVect");
342 }
343 
345  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
346 {
347  AssembleRHSElementVect(el, el, Tr, elvect);
348 }
349 
351  const FiniteElement &el1, const FiniteElement &el2,
352  FaceElementTransformations &Tr, Vector &elvect)
353 {
354  int dim, ndof1, ndof2, ndof, ndoftotal;
355  double w;
356  Vector temp_elvect;
357 
358  dim = el1.GetDim();
359  ndof1 = el1.GetDof();
360  ndof2 = el2.GetDof();
361  ndoftotal = ndof1 + ndof2;
362  if (Tr.Elem2No >= NEproc ||
364  {
365  ndoftotal = ndof1;
366  }
367 
368  elvect.SetSize(ndoftotal);
369  elvect = 0.0;
370 
371  bool elem1f = true;
372  int elem1 = Tr.Elem1No,
373  elem2 = Tr.Elem2No,
374  marker1 = (*elem_marker)[elem1];
375 
376  int marker2;
377  if (Tr.Elem2No >= NEproc)
378  {
379  marker2 = (*elem_marker)[NEproc+par_shared_face_count];
381  }
383  {
384  marker2 = marker1;
385  }
386  else
387  {
388  marker2 = (*elem_marker)[elem2];
389  }
390 
391  if (!include_cut_cell)
392  {
393  // 1 is inside and 2 is cut or 1 is a boundary element.
394  if ( marker1 == ShiftedFaceMarker::SBElementType::INSIDE &&
395  (marker2 == ls_cut_marker ||
397  {
398  elem1f = true;
399  ndof = ndof1;
400  }
401  // 1 is cut, 2 is inside
402  else if (marker1 == ls_cut_marker &&
403  marker2 == ShiftedFaceMarker::SBElementType::INSIDE)
404  {
405  if (Tr.Elem2No >= NEproc) { return; }
406  elem1f = false;
407  ndof = ndof2;
408  }
409  else
410  {
411  return;
412  }
413  }
414  else
415  {
416  // 1 is cut and 2 is outside or 1 is a boundary element.
417  if (marker1 == ls_cut_marker &&
418  (marker2 == ShiftedFaceMarker::SBElementType::OUTSIDE ||
420  {
421  elem1f = true;
422  ndof = ndof1;
423  }
424  // 1 is outside, 2 is cut
425  else if (marker1 == ShiftedFaceMarker::SBElementType::OUTSIDE &&
426  marker2 == ls_cut_marker)
427  {
428  if (Tr.Elem2No >= NEproc) { return; }
429  elem1f = false;
430  ndof = ndof2;
431  }
432  else
433  {
434  return;
435  }
436  }
437 
438  int offset = elem1f ? 0 : ndof1;
439  temp_elvect.SetDataAndSize(elvect.GetData()+offset, ndof);
440 
441  nor.SetSize(dim);
442  nh.SetSize(dim);
443  ni.SetSize(dim);
444  adjJ.SetSize(dim);
445 
446  shape.SetSize(ndof);
447  dshape.SetSize(ndof, dim);
448  dshape_dd.SetSize(ndof);
449  dshape_dn.SetSize(ndof);
450 
451  const IntegrationRule *ir = IntRule;
452  if (ir == NULL)
453  {
454  // a simple choice for the integration order; is this OK?
455  int order = elem1f ? 4*el1.GetOrder() : 4*el2.GetOrder();
456  ir = &IntRules.Get(Tr.GetGeometryType(), order);
457  }
458 
459  Array<DenseMatrix *> dkphi_dxk;
460  DenseMatrix grad_phys;
461  Vector Factorial;
462  Array<DenseMatrix *> grad_phys_dir;
463 
464  if (nterms > 0)
465  {
466  if (elem1f)
467  {
468  el1.ProjectGrad(el1, *Tr.Elem1, grad_phys);
469  }
470  else
471  {
472  el2.ProjectGrad(el2, *Tr.Elem2, grad_phys);
473  }
474 
475  DenseMatrix grad_work;
476  grad_phys_dir.SetSize(dim); // NxN matrix for derivative in each direction
477  for (int i = 0; i < dim; i++)
478  {
479  grad_phys_dir[i] = new DenseMatrix(ndof, ndof);
480  grad_phys_dir[i]->CopyRows(grad_phys, i*ndof, (i+1)*ndof-1);
481  }
482 
483 
484  DenseMatrix grad_phys_work = grad_phys;
485  grad_phys_work.SetSize(ndof, ndof*dim);
486 
487  dkphi_dxk.SetSize(nterms);
488 
489  for (int i = 0; i < nterms; i++)
490  {
491  int sz1 = pow(dim, i+1);
492  dkphi_dxk[i] = new DenseMatrix(ndof, ndof*sz1*dim);
493  int loc_col_per_dof = sz1;
494  for (int k = 0; k < dim; k++)
495  {
496  grad_work.SetSize(ndof, ndof*sz1);
497  // grad_work[k] has derivative in kth direction for each DOF.
498  // grad_work[0] has d^2phi/dx^2 and d^2phi/dxdy terms and
499  // grad_work[1] has d^2phi/dydx and d^2phi/dy2 terms for each dof
500  if (i == 0)
501  {
502  Mult(*grad_phys_dir[k], grad_phys_work, grad_work);
503  }
504  else
505  {
506  Mult(*grad_phys_dir[k], *dkphi_dxk[i-1], grad_work);
507  }
508  // Now we must place columns for each dof together so that they are
509  // in order: d^2phi/dx^2, d^2phi/dxdy, d^2phi/dydx, d^2phi/dy2.
510  for (int j = 0; j < ndof; j++)
511  {
512  for (int d = 0; d < loc_col_per_dof; d++)
513  {
514  Vector col;
515  int tot_col_per_dof = loc_col_per_dof*dim;
516  grad_work.GetColumn(j*loc_col_per_dof+d, col);
517  dkphi_dxk[i]->SetCol(j*tot_col_per_dof+k*loc_col_per_dof+d, col);
518  }
519  }
520  }
521  }
522 
523  for (int i = 0; i < grad_phys_dir.Size(); i++)
524  {
525  delete grad_phys_dir[i];
526  }
527 
528  Factorial.SetSize(nterms);
529  Factorial(0) = 2;
530  for (int i = 1; i < nterms; i++)
531  {
532  Factorial(i) = Factorial(i-1)*(i+2);
533  }
534  }
535 
536  DenseMatrix q_hess_dn(dim, ndof);
537  Vector q_hess_dn_work(q_hess_dn.GetData(), ndof*dim);
538  Vector q_hess_dot_d(ndof);
539 
540  Vector D(vD->GetVDim());
541  Vector wrk = shape;
542  // Assemble: -< u_D, \nabla w.n >
543  // +<alpha h^{-1} u_D, w + \nabla w.d + h.o.t>
544  for (int p = 0; p < ir->GetNPoints(); p++)
545  {
546  const IntegrationPoint &ip = ir->IntPoint(p);
547 
548  // Set the integration point in the face and the neighboring element
549  Tr.SetAllIntPoints(&ip);
550 
551  // Access the neighboring element's integration point
552  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
553  const IntegrationPoint &eip1 = Tr.GetElement1IntPoint();
554  const IntegrationPoint &eip2 = Tr.GetElement2IntPoint();
555 
556  if (dim == 1)
557  {
558  nor(0) = 2*eip.x - 1.0;
559  }
560  else
561  {
562  CalcOrtho(Tr.Jacobian(), nor);
563  }
564  vD->Eval(D, Tr, ip);
565 
566  // Make sure the normal vector is pointing outside the domain.
567  if (!elem1f) { nor *= -1; }
568 
569  double hinvdx;
570 
571  if (elem1f)
572  {
573  el1.CalcShape(eip1, shape);
574  el1.CalcDShape(eip1, dshape);
575  hinvdx =nor*nor/Tr.Elem1->Weight();
576  w = ip.weight * uD->Eval(Tr, ip, D) / Tr.Elem1->Weight();
578  }
579  else
580  {
581  el2.CalcShape(eip2, shape);
582  el2.CalcDShape(eip2, dshape);
583  hinvdx = nor*nor/Tr.Elem2->Weight();
584  w = ip.weight * uD->Eval(Tr, ip, D) / Tr.Elem2->Weight();
586  }
587 
588  ni.Set(w, nor);
589  adjJ.Mult(ni, nh);
590 
592  temp_elvect.Add(-1., dshape_dn); // T2
593 
594  double jinv;
595  if (elem1f)
596  {
597  w = ip.weight * uD->Eval(Tr, ip, D) * alpha * hinvdx;
598  jinv = 1./Tr.Elem1->Weight();
599  }
600  else
601  {
602  w = ip.weight * uD->Eval(Tr, ip, D) * alpha * hinvdx;
603  jinv = 1./Tr.Elem2->Weight();
604  }
605  adjJ.Mult(D, nh);
606  nh *= jinv;
608 
609  q_hess_dot_d = 0.;
610  for (int i = 0; i < nterms; i++)
611  {
612  int sz1 = pow(dim, i+1);
613  DenseMatrix T1(dim, ndof*sz1);
614  Vector T1_wrk(T1.GetData(), dim*ndof*sz1);
615  dkphi_dxk[i]->MultTranspose(shape, T1_wrk);
616 
617  DenseMatrix T2;
618  Vector T2_wrk;
619  for (int j = 0; j < i+1; j++)
620  {
621  int sz2 = pow(dim, i-j);
622  T2.SetSize(dim, ndof*sz2);
623  T2_wrk.SetDataAndSize(T2.GetData(), dim*ndof*sz2);
624  T1.MultTranspose(D, T2_wrk);
625  T1 = T2;
626  }
627  Vector q_hess_dot_d_work(ndof);
628  T1.MultTranspose(D, q_hess_dot_d_work);
629  q_hess_dot_d_work *= 1./Factorial(i);
630  q_hess_dot_d += q_hess_dot_d_work;
631  }
632 
633  wrk = shape;
634  wrk += dshape_dd; // \grad w .d
635  wrk += q_hess_dot_d;
636  temp_elvect.Add(w, wrk); // <u, gradw.d>
637  }
638 
639  for (int i = 0; i < dkphi_dxk.Size(); i++)
640  {
641  delete dkphi_dxk[i];
642  }
643 }
644 
645 
647  const FiniteElement &el1, const FiniteElement &el2,
649 {
650  int dim, ndof1, ndof2, ndof, ndoftotal;
651  double w;
652  DenseMatrix temp_elmat;
653 
654  dim = el1.GetDim();
655  ndof1 = el1.GetDof();
656  ndof2 = el2.GetDof();
657  ndoftotal = Trans.ElementType == ElementTransformation::BDR_FACE ?
658  ndof1 : ndof1 + ndof2;
659 
660  elmat.SetSize(ndoftotal);
661  elmat = 0.0;
662 
663  bool elem1f = true; // flag indicating whether Trans.Elem1No is part of the
664  // surrogate domain or not.
665  int elem1 = Trans.Elem1No,
666  elem2 = Trans.Elem2No,
667  marker1 = (*elem_marker)[elem1];
668 
669  int marker2;
670 
671  if (Trans.Elem2No >= NEproc)
672  {
673  marker2 = (*elem_marker)[NEproc+par_shared_face_count];
675  }
677  {
678  marker2 = marker1;
679  }
680  else
681  {
682  marker2 = (*elem_marker)[elem2];
683  }
684 
685  if (!include_cut_cell)
686  {
687  // 1 is inside and 2 is cut or 1 is a boundary element.
688  if (marker1 == ShiftedFaceMarker::SBElementType::INSIDE &&
689  (cut_marker.Find(marker2) != -1 ||
691  {
692  elem1f = true;
693  }
694  // 1 is cut, 2 is inside
695  else if (cut_marker.Find(marker1) != -1 &&
696  marker2 == ShiftedFaceMarker::SBElementType::INSIDE)
697  {
698  if (Trans.Elem2No >= NEproc) { return; }
699  elem1f = false;
700  }
701  else
702  {
703  return;
704  }
705  }
706  else
707  {
708  // 1 is cut and 2 is outside or 1 is a boundary element.
709  if (cut_marker.Find(marker1) != -1 &&
710  (marker2 == ShiftedFaceMarker::SBElementType::OUTSIDE ||
712  {
713  elem1f = true;
714  }
715  // 1 is outside, 2 is cut
716  else if (marker1 == ShiftedFaceMarker::SBElementType::OUTSIDE &&
717  cut_marker.Find(marker2) != -1)
718  {
719  if (Trans.Elem2No >= NEproc) { return; }
720  elem1f = false;
721  }
722  else
723  {
724  return;
725  }
726  }
727 
728  ndof = elem1f ? ndof1 : ndof2;
729 
730  temp_elmat.SetSize(ndof);
731  temp_elmat = 0.;
732 
733  nor.SetSize(dim);
734  nh.SetSize(dim);
735  ni.SetSize(dim);
736  adjJ.SetSize(dim);
737 
738  shape.SetSize(ndof);
739  dshape.SetSize(ndof, dim);
740  dshapedn.SetSize(ndof);
741  Vector wrk = shape;
742 
743  const IntegrationRule *ir = IntRule;
744  if (ir == NULL)
745  {
746  int order = elem1f ? 4*el1.GetOrder() : 4*el2.GetOrder();
747  ir = &IntRules.Get(Trans.GetGeometryType(), order);
748  }
749 
750  Array<DenseMatrix *> dkphi_dxk;
751  DenseMatrix grad_phys;
752  Vector Factorial;
753  Array<DenseMatrix *> grad_phys_dir;
754 
755  if (nterms > 0)
756  {
757  if (elem1f)
758  {
759  el1.ProjectGrad(el1, *Trans.Elem1, grad_phys);
760  }
761  else
762  {
763  el2.ProjectGrad(el2, *Trans.Elem2, grad_phys);
764  }
765 
766  DenseMatrix grad_work;
767  grad_phys_dir.SetSize(dim); // NxN matrices for derivative in each direction
768  for (int i = 0; i < dim; i++)
769  {
770  grad_phys_dir[i] = new DenseMatrix(ndof, ndof);
771  grad_phys_dir[i]->CopyRows(grad_phys, i*ndof, (i+1)*ndof-1);
772  }
773 
774  DenseMatrix grad_phys_work = grad_phys;
775  grad_phys_work.SetSize(ndof, ndof*dim);
776 
777  dkphi_dxk.SetSize(nterms);
778 
779  for (int i = 0; i < nterms; i++)
780  {
781  int sz1 = pow(dim, i+1);
782  dkphi_dxk[i] = new DenseMatrix(ndof, ndof*sz1*dim);
783  int loc_col_per_dof = sz1;
784  int tot_col_per_dof = loc_col_per_dof*dim;
785  for (int k = 0; k < dim; k++)
786  {
787  grad_work.SetSize(ndof, ndof*sz1);
788  // grad_work[k] has derivative in kth direction for each DOF.
789  // grad_work[0] has d^2phi/dx^2 and d^2phi/dxdy terms and
790  // grad_work[1] has d^2phi/dydx and d^2phi/dy2 terms for each dof
791  if (i == 0)
792  {
793  Mult(*grad_phys_dir[k], grad_phys_work, grad_work);
794  }
795  else
796  {
797  Mult(*grad_phys_dir[k], *dkphi_dxk[i-1], grad_work);
798  }
799  // Now we must place columns for each dof together so that they are
800  // in order: d^2phi/dx^2, d^2phi/dxdy, d^2phi/dydx, d^2phi/dy2.
801  for (int j = 0; j < ndof; j++)
802  {
803  for (int d = 0; d < loc_col_per_dof; d++)
804  {
805  Vector col;
806  grad_work.GetColumn(j*loc_col_per_dof+d, col);
807  dkphi_dxk[i]->SetCol(j*tot_col_per_dof+k*loc_col_per_dof+d, col);
808  }
809  }
810  }
811  }
812 
813  for (int i = 0; i < grad_phys_dir.Size(); i++)
814  {
815  delete grad_phys_dir[i];
816  }
817 
818  Factorial.SetSize(nterms);
819  Factorial(0) = 1;
820  for (int i = 1; i < nterms; i++)
821  {
822  Factorial(i) = Factorial(i-1)*(i+1);
823  }
824  }
825 
826 
827  DenseMatrix q_hess_dn(dim, ndof);
828  Vector q_hess_dn_work(q_hess_dn.GetData(), ndof*dim);
829  Vector q_hess_dot_d_nhat(ndof);
830 
831  Vector D(vD->GetVDim());
832  Vector Nhat(vN->GetVDim());
833  // Assemble: <nabla(nabla u).d.nhat (n.nhat), w>
834  for (int p = 0; p < ir->GetNPoints(); p++)
835  {
836  const IntegrationPoint &ip = ir->IntPoint(p);
837 
838  // Set the integration point in the face and the neighboring elements
839  Trans.SetAllIntPoints(&ip);
840 
841  // Access the neighboring elements' integration points
842  // Note: eip2 will only contain valid data if Elem2 exists
843  const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
844  const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
845 
846  if (dim == 1)
847  {
848  nor(0) = 2*eip1.x - 1.0;
849  }
850  else
851  {
852  // Note: this normal accounts for the weight of the surface transformation
853  // Jacobian i.e. nor = nhat*det(J)
854  CalcOrtho(Trans.Jacobian(), nor);
855  }
856  vD->Eval(D, Trans, ip);
857  vN->Eval(Nhat, Trans, ip, D);
858 
859  // Make sure the normal vector is pointing outside the domain.
860  if (!elem1f) { nor *= -1; }
861 
862  if (elem1f)
863  {
864  el1.CalcShape(eip1, shape);
865  el1.CalcDShape(eip1, dshape);
866  w = ip.weight/Trans.Elem1->Weight();
867  CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
868  }
869  else
870  {
871  el1.CalcShape(eip2, shape);
872  el1.CalcDShape(eip2, dshape);
873  w = ip.weight/Trans.Elem2->Weight();
874  CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
875  }
876 
877  ni.Set(w, nor); // nor/det(J)
878  adjJ.Mult(ni, nh);
879  dshape.Mult(nh, dshapedn); //dphi/dn * Jinv * nor
880 
881  // -<w, grad u.n> - Term 2
882  AddMult_a_VWt(-1., shape, dshapedn, temp_elmat);
883 
884  // <w, (grad u.nhat)nhat.n
885  // Here nhat is the normal vector at true boundary and n is the normal
886  // vector at shifted boundary
887  double n_dot_ntilde = nor*Nhat;
888  ni.Set(w, Nhat);
889  adjJ.Mult(ni, nh);
891  dshapedn *= n_dot_ntilde;
892  AddMult_a_VWt(1., shape, dshapedn, temp_elmat);
893 
894  q_hess_dot_d_nhat = 0.;
895  for (int i = 0; i < nterms; i++)
896  {
897  int sz1 = pow(dim, i+1);
898  DenseMatrix T1(dim, ndof*sz1);
899  Vector T1_wrk(T1.GetData(), dim*ndof*sz1);
900  dkphi_dxk[i]->MultTranspose(shape, T1_wrk);
901 
902  DenseMatrix T2;
903  Vector T2_wrk;
904  for (int j = 0; j < i+1; j++)
905  {
906  int sz2 = pow(dim, i-j);
907  T2.SetSize(dim, ndof*sz2);
908  T2_wrk.SetDataAndSize(T2.GetData(), dim*ndof*sz2);
909  T1.MultTranspose(D, T2_wrk);
910  T1 = T2;
911  }
912  Vector q_hess_dot_d_work(ndof);
913  T1.MultTranspose(Nhat, q_hess_dot_d_work);
914  q_hess_dot_d_work *= 1./Factorial(i);
915  q_hess_dot_d_nhat += q_hess_dot_d_work;
916  }
917 
918  wrk = q_hess_dot_d_nhat;
919  wrk *= ip.weight * n_dot_ntilde;
920 
921  AddMult_a_VWt(1., shape, wrk, temp_elmat);
922  } //p < ir->GetNPoints()
923  int offset = elem1f ? 0 : ndof1;
924  elmat.CopyMN(temp_elmat, offset, offset);
925 
926  for (int i = 0; i < dkphi_dxk.Size(); i++)
927  {
928  delete dkphi_dxk[i];
929  }
930 }
931 
933  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
934 {
935  mfem_error("SBM2NeumannLFIntegrator::AssembleRHSElementVect");
936 }
937 
939  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
940 {
941  AssembleRHSElementVect(el, el, Tr, elvect);
942 }
943 
945  const FiniteElement &el1, const FiniteElement &el2,
946  FaceElementTransformations &Tr, Vector &elvect)
947 {
948  int dim, ndof1, ndof2, ndof, ndoftotal;
949  double w;
950  Vector temp_elvect;
951 
952  dim = el1.GetDim();
953  ndof1 = el1.GetDof();
954  ndof2 = el2.GetDof();
955  ndoftotal = ndof1 + ndof2;
956  if (Tr.Elem2No >= NEproc ||
958  {
959  ndoftotal = ndof1;
960  }
961 
962  elvect.SetSize(ndoftotal);
963  elvect = 0.0;
964 
965  bool elem1f = true;
966  int elem1 = Tr.Elem1No,
967  elem2 = Tr.Elem2No,
968  marker1 = (*elem_marker)[elem1];
969 
970  int marker2;
971  if (Tr.Elem2No >= NEproc)
972  {
973  marker2 = (*elem_marker)[NEproc+par_shared_face_count];
975  }
977  {
978  marker2 = marker1;
979  }
980  else
981  {
982  marker2 = (*elem_marker)[elem2];
983  }
984  if (!include_cut_cell)
985  {
986  // 1 is inside and 2 is cut or 1 is a boundary element.
987  if ( marker1 == ShiftedFaceMarker::SBElementType::INSIDE &&
988  (marker2 == ls_cut_marker ||
990  {
991  elem1f = true;
992  ndof = ndof1;
993  }
994  // 1 is cut, 2 is inside
995  else if (marker1 == ls_cut_marker &&
996  marker2 == ShiftedFaceMarker::SBElementType::INSIDE)
997  {
998  if (Tr.Elem2No >= NEproc) { return; }
999  elem1f = false;
1000  ndof = ndof2;
1001  }
1002  else
1003  {
1004  return;
1005  }
1006  }
1007  else
1008  {
1009  // 1 is cut and 2 is outside or 1 is a boundary element.
1010  if (marker1 == ls_cut_marker &&
1011  (marker2 == ShiftedFaceMarker::SBElementType::OUTSIDE ||
1013  {
1014  elem1f = true;
1015  ndof = ndof1;
1016  }
1017  // 1 is outside, 2 is cut
1018  else if (marker1 == ShiftedFaceMarker::SBElementType::OUTSIDE &&
1019  marker2 == ls_cut_marker)
1020  {
1021  if (Tr.Elem2No >= NEproc) { return; }
1022  elem1f = false;
1023  ndof = ndof2;
1024  }
1025  else
1026  {
1027  return;
1028  }
1029  }
1030 
1031  int offset = elem1f ? 0 : ndof1;
1032  temp_elvect.SetDataAndSize(elvect.GetData()+offset, ndof);
1033 
1034  nor.SetSize(dim);
1035  shape.SetSize(ndof);
1036 
1037  const IntegrationRule *ir = IntRule;
1038  if (ir == NULL)
1039  {
1040  // a simple choice for the integration order; is this OK?
1041  int order = elem1f ? 4*el1.GetOrder() : 4*el2.GetOrder();
1042  ir = &IntRules.Get(Tr.GetGeometryType(), order);
1043  }
1044 
1045  Vector D(vD->GetVDim());
1046  Vector Nhat(vN->GetVDim());
1047  Vector wrk = shape;
1048  for (int p = 0; p < ir->GetNPoints(); p++)
1049  {
1050  const IntegrationPoint &ip = ir->IntPoint(p);
1051 
1052  // Set the integration point in the face and the neighboring element
1053  Tr.SetAllIntPoints(&ip);
1054 
1055  // Access the neighboring element's integration point
1056  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
1057  const IntegrationPoint &eip1 = Tr.GetElement1IntPoint();
1058  const IntegrationPoint &eip2 = Tr.GetElement2IntPoint();
1059 
1060  if (dim == 1)
1061  {
1062  nor(0) = 2*eip.x - 1.0;
1063  }
1064  else
1065  {
1066  CalcOrtho(Tr.Jacobian(), nor);
1067  }
1068  vD->Eval(D, Tr, ip);
1069  vN->Eval(Nhat, Tr, ip, D);
1070 
1071  // Make sure the normal vector is pointing outside the domain.
1072  if (!elem1f) { nor *= -1; }
1073 
1074  if (elem1f)
1075  {
1076  el1.CalcShape(eip1, shape);
1077  w = ip.weight * uN->Eval(Tr, ip, D);
1078  }
1079  else
1080  {
1081  el2.CalcShape(eip2, shape);
1082  w = ip.weight * uN->Eval(Tr, ip, D);
1083  }
1084 
1085  double n_dot_ntilde = nor*Nhat;
1086  wrk.Set(n_dot_ntilde*w, shape);
1087  //<w, (nhat.n)t_n)
1088  temp_elvect.Add(1., wrk);
1089  }
1090 }
1091 
1092 }
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
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:162
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
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:472
ShiftedVectorFunctionCoefficient * vN
Definition: sbm_solver.hpp:273
VectorCoefficient * vD
Definition: sbm_solver.hpp:92
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:2074
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
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
const IntegrationPoint & GetElement2IntPoint()
Get a const reference to the integration point in neighboring element 2 corresponding to the currentl...
Definition: eltrans.hpp:587
VectorCoefficient * vD
Definition: sbm_solver.hpp:275
const IntegrationRule * IntRule
Definition: lininteg.hpp:25
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
ElementTransformation * Elem2
Definition: eltrans.hpp:522
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:160
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2288
std::function< void(const Vector &, Vector &)> Function
Definition: sbm_solver.hpp:54
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:188
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
ShiftedFunctionCoefficient * uN
Definition: sbm_solver.hpp:274
void AddMult_a_VWt(const double a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
Definition: densemat.cpp:2806
void SetAllIntPoints(const IntegrationPoint *face_ip)
Set the integration point in the Face and the two neighboring elements, if present.
Definition: eltrans.hpp:566
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:119
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: sbm_solver.cpp:338
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
virtual void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
Definition: sbm_solver.cpp:646
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
int GetVDim()
Returns dimension of the vector.
ShiftedVectorFunctionCoefficient * vN
Definition: sbm_solver.hpp:211
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Definition: sbm_solver.hpp:62
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1270
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
VectorCoefficient * vD
Definition: sbm_solver.hpp:212
int Find(const T &el) const
Return the first index where &#39;el&#39; is found; return -1 if not found.
Definition: array.hpp:810
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
const IntegrationPoint & GetElement1IntPoint()
Get a const reference to the integration point in neighboring element 1 corresponding to the currentl...
Definition: eltrans.hpp:577
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
Definition: sbm_solver.hpp:36
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:156
const IntegrationRule * IntRule
Definition: nonlininteg.hpp:30
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:267
std::function< double(const Vector &)> Function
Definition: sbm_solver.hpp:26
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
int dim
Definition: ex24.cpp:53
virtual void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
Definition: sbm_solver.cpp:48
ElementTransformation * Elem1
Definition: eltrans.hpp:522
Vector data type.
Definition: vector.hpp:60
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:160
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco)
Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this.
Definition: densemat.cpp:1517
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: sbm_solver.cpp:932
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
ShiftedFunctionCoefficient * uD
Definition: sbm_solver.hpp:151