MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
paramnonlinearform.cpp
Go to the documentation of this file.
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 
13 #include "mfem.hpp"
14 #include "paramnonlinearform.hpp"
15 
16 
17 namespace mfem
18 {
19 
20 double ParametricBNLFormIntegrator::GetElementEnergy(const
21  Array<const FiniteElement *> &el,
22  const Array<const FiniteElement *> &pel,
23  ElementTransformation &Tr,
24  const Array<const Vector *> &elfun,
25  const Array<const Vector *> &pelfun)
26 {
27  mfem_error("ParametricBNLFormIntegrator::GetElementEnergy"
28  " is not overloaded!");
29  return 0.0;
30 }
31 
32 void ParametricBNLFormIntegrator::AssembleFaceGrad(const
33  Array<const FiniteElement *> &el1,
34  const Array<const FiniteElement *> &el2,
35  const Array<const FiniteElement *> &pel1,
36  const Array<const FiniteElement *> &pel2,
37  FaceElementTransformations &Tr,
38  const Array<const Vector *> &elfun,
39  const Array<const Vector *> &pelfun,
40  const Array2D<DenseMatrix *> &elmats)
41 {
42  mfem_error("ParametricBNLFormIntegrator::AssembleFaceGrad"
43  " is not overloaded!");
44 }
45 
46 void ParametricBNLFormIntegrator::AssembleElementGrad(const
47  Array<const FiniteElement *> &el,
48  const Array<const FiniteElement *> &pel,
49  ElementTransformation &Tr,
50  const Array<const Vector *> &elfun,
51  const Array<const Vector *> &pelfun,
52  const Array2D<DenseMatrix *> &elmats)
53 {
54  mfem_error("ParametricBNLFormIntegrator::AssembleElementGrad"
55  " is not overloaded!");
56 }
57 
58 void ParametricBNLFormIntegrator::AssembleElementVector(
59  const Array<const FiniteElement *> &el,
60  const Array<const FiniteElement *> &pel,
61  ElementTransformation &Tr,
62  const Array<const Vector *> &elfun,
63  const Array<const Vector *> &pelfun,
64  const Array<Vector *> &elvec)
65 {
66  mfem_error("ParametricBNLFormIntegrator::AssembleElementVector"
67  " is not overloaded!");
68 }
69 
70 void ParametricBNLFormIntegrator::AssembleFaceVector(const
71  Array<const FiniteElement *> &el1,
72  const Array<const FiniteElement *> &el2,
73  const Array<const FiniteElement *> &pel1,
74  const Array<const FiniteElement *> &pel2,
75  FaceElementTransformations &Tr,
76  const Array<const Vector *> &elfun,
77  const Array<const Vector *> &pelfun,
78  const Array<Vector *> &elvect)
79 {
80  mfem_error("ParametricBNLFormIntegrator::AssembleFaceVector"
81  " is not overloaded!");
82 }
83 
84 
85 void ParametricBNLFormIntegrator::AssemblePrmElementVector(
86  const Array<const FiniteElement *> &el,
87  const Array<const FiniteElement *> &pel,
88  ElementTransformation &Tr,
89  const Array<const Vector *> &elfun,
90  const Array<const Vector *> &alfun,
91  const Array<const Vector *> &pelfun,
92  const Array<Vector *> &elvec)
93 {
94  mfem_error("ParametricBNLFormIntegrator::AssemblePrmElementVector"
95  " is not overloaded!");
96 }
97 
98 void ParametricBNLFormIntegrator::AssemblePrmFaceVector(
99  const Array<const FiniteElement *> &el1,
100  const Array<const FiniteElement *> &el2,
101  const Array<const FiniteElement *> &pel1,
102  const Array<const FiniteElement *> &pel2,
103  FaceElementTransformations &Tr,
104  const Array<const Vector *> &elfun,
105  const Array<const Vector *> &alfun,
106  const Array<const Vector *> &pelfun,
107  const Array<Vector *> &elvect)
108 {
109  mfem_error("ParametricBNLFormIntegrator::AssemblePrmFaceVector"
110  " is not overloaded!");
111 }
112 
113 
114 ParametricBNLForm::ParametricBNLForm(): fes(0), paramfes(0), BlockGrad(nullptr)
115 {
116  height = 0;
117  width = 0;
118 
119  paramheight = 0;
120  paramwidth = 0;
121 }
122 
123 ParametricBNLForm::ParametricBNLForm(Array<FiniteElementSpace *> &statef,
124  Array<FiniteElementSpace *> &paramf):
125  fes(0),paramfes(0), BlockGrad(nullptr)
126 {
127  SetSpaces(statef,paramf);
128 }
129 
130 void ParametricBNLForm::SetSpaces(Array<FiniteElementSpace *> &statef,
131  Array<FiniteElementSpace *> &paramf)
132 {
133  delete BlockGrad;
134  BlockGrad = nullptr;
135  for (int i=0; i<Grads.NumRows(); ++i)
136  {
137  for (int j=0; j<Grads.NumCols(); ++j)
138  {
139  delete Grads(i,j);
140  delete cGrads(i,j);
141  }
142  }
143  for (int i = 0; i < ess_tdofs.Size(); ++i)
144  {
145  delete ess_tdofs[i];
146  }
147  for (int i = 0; i < paramess_tdofs.Size(); ++i)
148  {
149  delete paramess_tdofs[i];
150  }
151 
152  height = 0;
153  width = 0;
154  // set state feses
155  statef.Copy(fes);
156  // set the block sizes
157  block_offsets.SetSize(statef.Size() + 1);
158  block_trueOffsets.SetSize(statef.Size() + 1);
159  block_offsets[0] = 0;
160  block_trueOffsets[0] = 0;
161  for (int i=0; i<fes.Size(); ++i)
162  {
163  block_offsets[i+1] = fes[i]->GetVSize();
164  block_trueOffsets[i+1] = fes[i]->GetTrueVSize();
165  }
166  block_offsets.PartialSum();
167  block_trueOffsets.PartialSum();
168 
169  // Set design/parametric feses
170  paramf.Copy(paramfes);
171  paramblock_offsets.SetSize(paramf.Size() + 1);
172  paramblock_trueOffsets.SetSize(paramf.Size() + 1);
173  paramblock_offsets[0] = 0;
174  paramblock_trueOffsets[0] = 0;
175  for (int i=0; i<paramfes.Size(); ++i)
176  {
177  paramblock_offsets[i+1] = paramfes[i]->GetVSize();
178  paramblock_trueOffsets[i+1] = paramfes[i]->GetTrueVSize();
179  }
180  paramblock_offsets.PartialSum();
181  paramblock_trueOffsets.PartialSum();
182 
183  // Set the size of the operator
184  height = block_trueOffsets[fes.Size()];
185  width = block_trueOffsets[fes.Size()];
186 
187  Grads.SetSize(fes.Size(), fes.Size());
188  Grads = NULL;
189 
190  cGrads.SetSize(fes.Size(), fes.Size());
191  cGrads = NULL;
192 
193  P.SetSize(fes.Size());
194  cP.SetSize(fes.Size());
195  ess_tdofs.SetSize(fes.Size());
196 
197  for (int s = 0; s < fes.Size(); ++s)
198  {
199  // Retrieve prolongation matrix for each FE space
200  P[s] = fes[s]->GetProlongationMatrix();
201  cP[s] = dynamic_cast<const SparseMatrix *>(P[s]);
202 
203  // If the P Operator exists and its type is not SparseMatrix, this
204  // indicates the Operator is part of parallel run.
205  if (P[s] && !cP[s])
206  {
207  is_serial = false;
208  }
209 
210  // If the P Operator exists and its type is SparseMatrix, this indicates
211  // the Operator is serial but needs prolongation on assembly.
212  if (cP[s])
213  {
214  needs_prolongation = true;
215  }
216 
217  ess_tdofs[s] = new Array<int>;
218  }
219 
220  // Set the size of the design/parametric part.
221  paramheight = paramblock_trueOffsets[paramfes.Size()];
222  paramwidth = paramblock_trueOffsets[paramfes.Size()];
223 
224  Pparam.SetSize(paramfes.Size());
225  cPparam.SetSize(paramfes.Size());
226  paramess_tdofs.SetSize(paramfes.Size());
227 
228  for (int s = 0; s < paramfes.Size(); ++s)
229  {
230  // Retrieve prolongation matrix for each FE space.
231  Pparam[s] = paramfes[s]->GetProlongationMatrix();
232  cPparam[s] = dynamic_cast<const SparseMatrix *>(Pparam[s]);
233 
234  // If the P Operator exists and its type is SparseMatrix, this indicates
235  // the Operator is serial but needs prolongation on assembly.
236  if (cPparam[s])
237  {
238  prmneeds_prolongation = true;
239  }
240 
241  paramess_tdofs[s] = new Array<int>;
242  }
243 
244  xsv.Update(block_offsets);
245  adv.Update(block_offsets);
246  xdv.Update(paramblock_offsets);
247 
248 }
249 
250 void ParametricBNLForm::AddBdrFaceIntegrator(ParametricBNLFormIntegrator *nlfi,
251  Array<int> &bdr_marker)
252 {
253  bfnfi.Append(nlfi);
254  bfnfi_marker.Append(&bdr_marker);
255 }
256 
257 double ParametricBNLForm::GetEnergyBlocked(const BlockVector &bx,
258  const BlockVector &dx) const
259 {
260  Array<Array<int> *> vdofs(fes.Size());
261  Array<Vector *> el_x(fes.Size());
262  Array<const Vector *> el_x_const(fes.Size());
263  Array<const FiniteElement *> fe(fes.Size());
264  ElementTransformation *T;
265 
266  Array<Array<int> *> prmvdofs(paramfes.Size());
267  Array<Vector *> prmel_x(paramfes.Size());
268  Array<const Vector *> prmel_x_const(paramfes.Size());
269  Array<const FiniteElement *> prmfe(paramfes.Size());
270 
271  double energy = 0.0;
272 
273  for (int i=0; i<fes.Size(); ++i)
274  {
275  el_x_const[i] = el_x[i] = new Vector();
276  vdofs[i] = new Array<int>;
277  }
278 
279  for (int i=0; i<paramfes.Size(); ++i)
280  {
281  prmel_x_const[i] = prmel_x[i] = new Vector();
282  prmvdofs[i] = new Array<int>;
283  }
284 
285 
286  if (dnfi.Size())
287  {
288  for (int i = 0; i < fes[0]->GetNE(); ++i)
289  {
290  T = fes[0]->GetElementTransformation(i);
291  for (int s=0; s<fes.Size(); ++s)
292  {
293  fe[s] = fes[s]->GetFE(i);
294  fes[s]->GetElementVDofs(i, *vdofs[s]);
295  bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
296  }
297 
298  for (int s=0; s<paramfes.Size(); ++s)
299  {
300  prmfe[s] = paramfes[s]->GetFE(i);
301  paramfes[s]->GetElementVDofs(i, *prmvdofs[s]);
302  dx.GetBlock(s).GetSubVector(*prmvdofs[s], *prmel_x[s]);
303  }
304 
305 
306 
307  for (int k = 0; k < dnfi.Size(); ++k)
308  {
309  energy += dnfi[k]->GetElementEnergy(fe, prmfe, *T, el_x_const, prmel_x_const);
310  }
311  }
312  }
313 
314  // Free the allocated memory
315  for (int i = 0; i < fes.Size(); ++i)
316  {
317  delete el_x[i];
318  delete vdofs[i];
319  }
320 
321  for (int i = 0; i < paramfes.Size(); ++i)
322  {
323  delete prmel_x[i];
324  delete prmvdofs[i];
325  }
326 
327  if (fnfi.Size())
328  {
329  MFEM_ABORT("TODO: add energy contribution from interior face terms");
330  }
331 
332  if (bfnfi.Size())
333  {
334  MFEM_ABORT("TODO: add energy contribution from boundary face terms");
335  }
336 
337  return energy;
338 }
339 
340 void ParametricBNLForm::SetStateFields(const Vector &xv) const
341 {
342  BlockVector bx(const_cast<Vector&>(xv), block_trueOffsets);
343  if (needs_prolongation)
344  {
345  for (int s = 0; s < fes.Size(); s++)
346  {
347  P[s]->Mult(bx.GetBlock(s), xsv.GetBlock(s));
348  }
349  }
350  else
351  {
352  xsv=bx;
353  }
354 }
355 
356 void ParametricBNLForm::SetAdjointFields(const Vector &av) const
357 {
358  BlockVector bx(const_cast<Vector&>(av), block_trueOffsets);
359  if (needs_prolongation)
360  {
361  for (int s = 0; s < fes.Size(); s++)
362  {
363  P[s]->Mult(bx.GetBlock(s), adv.GetBlock(s));
364  }
365  }
366  else
367  {
368  adv=bx;
369  }
370 
371 }
372 
373 void ParametricBNLForm::SetParamFields(const Vector &dv) const
374 {
375  BlockVector bx(const_cast<Vector&>(dv), paramblock_trueOffsets);
376  if (prmneeds_prolongation)
377  {
378  for (int s = 0; s < paramfes.Size(); s++)
379  {
380  Pparam[s]->Mult(bx.GetBlock(s), xdv.GetBlock(s));
381  }
382  }
383  else
384  {
385  xdv=bx;
386  }
387 }
388 
389 double ParametricBNLForm::GetEnergy(const Vector &x) const
390 {
391  xs.Update(const_cast<Vector&>(x),block_offsets);
392  return GetEnergyBlocked(xs,xdv);
393 }
394 
395 void ParametricBNLForm::SetEssentialBC(const Array<Array<int> *>
396  &bdr_attr_is_ess,
397  Array<Vector *> &rhs)
398 {
399  for (int s = 0; s < fes.Size(); ++s)
400  {
401  ess_tdofs[s]->SetSize(ess_tdofs.Size());
402 
403  fes[s]->GetEssentialTrueDofs(*bdr_attr_is_ess[s], *ess_tdofs[s]);
404 
405  if (rhs[s])
406  {
407  rhs[s]->SetSubVector(*ess_tdofs[s], 0.0);
408  }
409  }
410 }
411 
412 void ParametricBNLForm::SetParamEssentialBC(const Array<Array<int> *>
413  &bdr_attr_is_ess,
414  Array<Vector *> &rhs)
415 {
416  for (int s = 0; s < paramfes.Size(); ++s)
417  {
418  paramess_tdofs[s]->SetSize(paramess_tdofs.Size());
419 
420  paramfes[s]->GetEssentialTrueDofs(*bdr_attr_is_ess[s], *paramess_tdofs[s]);
421 
422  if (rhs[s])
423  {
424  rhs[s]->SetSubVector(*paramess_tdofs[s], 0.0);
425  }
426  }
427 }
428 
429 void ParametricBNLForm::MultParamBlocked(const BlockVector &bx,
430  const BlockVector &ax,
431  const BlockVector &dx,
432  BlockVector &dy) const
433 {
434  // State fields
435  Array<Array<int> *>vdofs(fes.Size());
436  Array<Array<int> *>vdofs2(fes.Size());
437  Array<Vector *> el_x(fes.Size());
438  Array<const Vector *> el_x_const(fes.Size());
439  Array<const FiniteElement *> fe(fes.Size());
440  Array<const FiniteElement *> fe2(fes.Size());
441 
442  ElementTransformation *T;
443 
444  // Adjoint fields
445  Array<Vector *> ael_x(fes.Size());
446  Array<const Vector *> ael_x_const(fes.Size());
447 
448  // Parametric fields
449  Array<Array<int> *>prmvdofs(paramfes.Size());
450  Array<Array<int> *>prmvdofs2(paramfes.Size());
451  Array<Vector *> prmel_x(paramfes.Size());
452  Array<const Vector *> prmel_x_const(paramfes.Size());
453  Array<Vector *> prmel_y(paramfes.Size());
454  Array<const FiniteElement *> prmfe(paramfes.Size());
455  Array<const FiniteElement *> prmfe2(paramfes.Size());
456 
457  dy = 0.0;
458 
459  for (int s=0; s<fes.Size(); ++s)
460  {
461  el_x_const[s] = el_x[s] = new Vector();
462  vdofs[s] = new Array<int>;
463  vdofs2[s] = new Array<int>;
464 
465  ael_x_const[s] = ael_x[s] = new Vector();
466  }
467 
468  for (int s=0; s<paramfes.Size(); ++s)
469  {
470  prmel_x_const[s] = prmel_x[s] = new Vector();
471  prmel_y[s] = new Vector();
472  prmvdofs[s] = new Array<int>;
473  prmvdofs2[s] = new Array<int>;
474  }
475 
476  if (dnfi.Size())
477  {
478  for (int i = 0; i < fes[0]->GetNE(); ++i)
479  {
480  T = fes[0]->GetElementTransformation(i);
481  for (int s = 0; s < fes.Size(); ++s)
482  {
483  fes[s]->GetElementVDofs(i, *(vdofs[s]));
484  fe[s] = fes[s]->GetFE(i);
485  bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
486  ax.GetBlock(s).GetSubVector(*(vdofs[s]), *ael_x[s]);
487  }
488 
489  for (int s = 0; s < paramfes.Size(); ++s)
490  {
491  paramfes[s]->GetElementVDofs(i, *(prmvdofs[s]));
492  prmfe[s] = paramfes[s]->GetFE(i);
493  dx.GetBlock(s).GetSubVector(*(prmvdofs[s]), *prmel_x[s]);
494  }
495 
496 
497  for (int k = 0; k < dnfi.Size(); ++k)
498  {
499  dnfi[k]->AssemblePrmElementVector(fe,prmfe, *T,
500  el_x_const, ael_x_const, prmel_x_const,
501  prmel_y);
502 
503  for (int s=0; s<paramfes.Size(); ++s)
504  {
505  if (prmel_y[s]->Size() == 0) { continue; }
506  dy.GetBlock(s).AddElementVector(*(prmvdofs[s]), *prmel_y[s]);
507  }
508  }
509  }
510  }
511 
512  if (fnfi.Size())
513  {
514  Mesh *mesh = fes[0]->GetMesh();
515  FaceElementTransformations *tr;
516 
517  for (int i = 0; i < mesh->GetNumFaces(); ++i)
518  {
519  tr = mesh->GetInteriorFaceTransformations(i);
520  if (tr != NULL)
521  {
522  for (int s=0; s<fes.Size(); ++s)
523  {
524  fe[s] = fes[s]->GetFE(tr->Elem1No);
525  fe2[s] = fes[s]->GetFE(tr->Elem2No);
526 
527  fes[s]->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
528  fes[s]->GetElementVDofs(tr->Elem2No, *(vdofs2[s]));
529 
530  vdofs[s]->Append(*(vdofs2[s]));
531 
532  bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
533  ax.GetBlock(s).GetSubVector(*(vdofs[s]), *ael_x[s]);
534  }
535 
536  for (int s=0; s<paramfes.Size(); ++s)
537  {
538  prmfe[s] = paramfes[s]->GetFE(tr->Elem1No);
539  prmfe2[s] = paramfes[s]->GetFE(tr->Elem2No);
540 
541  paramfes[s]->GetElementVDofs(tr->Elem1No, *(prmvdofs[s]));
542  paramfes[s]->GetElementVDofs(tr->Elem2No, *(prmvdofs2[s]));
543 
544  prmvdofs[s]->Append(*(prmvdofs2[s]));
545 
546  dx.GetBlock(s).GetSubVector(*(prmvdofs[s]), *prmel_x[s]);
547  }
548 
549  for (int k = 0; k < fnfi.Size(); ++k)
550  {
551 
552  fnfi[k]->AssemblePrmFaceVector(fe, fe2, prmfe, prmfe2, *tr,
553  el_x_const, ael_x_const, prmel_x_const, prmel_y);
554 
555  for (int s=0; s<paramfes.Size(); ++s)
556  {
557  if (prmel_y[s]->Size() == 0) { continue; }
558  dy.GetBlock(s).AddElementVector(*(prmvdofs[s]), *prmel_y[s]);
559  }
560  }
561  }
562  }
563  }
564 
565  if (bfnfi.Size())
566  {
567  Mesh *mesh = fes[0]->GetMesh();
568  FaceElementTransformations *tr;
569  // Which boundary attributes need to be processed?
570  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
571  mesh->bdr_attributes.Max() : 0);
572  bdr_attr_marker = 0;
573  for (int k = 0; k < bfnfi.Size(); ++k)
574  {
575  if (bfnfi_marker[k] == NULL)
576  {
577  bdr_attr_marker = 1;
578  break;
579  }
580  Array<int> &bdr_marker = *bfnfi_marker[k];
581  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
582  "invalid boundary marker for boundary face integrator #"
583  << k << ", counting from zero");
584  for (int i = 0; i < bdr_attr_marker.Size(); ++i)
585  {
586  bdr_attr_marker[i] |= bdr_marker[i];
587  }
588  }
589 
590  for (int i = 0; i < mesh->GetNBE(); ++i)
591  {
592  const int bdr_attr = mesh->GetBdrAttribute(i);
593  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
594 
595  tr = mesh->GetBdrFaceTransformations(i);
596  if (tr != NULL)
597  {
598  for (int s=0; s<fes.Size(); ++s)
599  {
600  fe[s] = fes[s]->GetFE(tr->Elem1No);
601  fe2[s] = fes[s]->GetFE(tr->Elem1No);
602 
603  fes[s]->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
604  bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
605  ax.GetBlock(s).GetSubVector(*(vdofs[s]), *ael_x[s]);
606  }
607 
608  for (int s=0; s<paramfes.Size(); ++s)
609  {
610  prmfe[s] = paramfes[s]->GetFE(tr->Elem1No);
611  prmfe2[s] = paramfes[s]->GetFE(tr->Elem1No);
612 
613  paramfes[s]->GetElementVDofs(tr->Elem1No, *(prmvdofs[s]));
614  dx.GetBlock(s).GetSubVector(*(prmvdofs[s]), *prmel_x[s]);
615  }
616 
617 
618  for (int k = 0; k < bfnfi.Size(); ++k)
619  {
620  if (bfnfi_marker[k] &&
621  (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
622 
623  bfnfi[k]->AssemblePrmFaceVector(fe, fe2, prmfe, prmfe2, *tr,
624  el_x_const, ael_x_const, prmel_x_const, prmel_y);
625 
626  for (int s=0; s<paramfes.Size(); ++s)
627  {
628  if (prmel_y[s]->Size() == 0) { continue; }
629  dy.GetBlock(s).AddElementVector(*(prmvdofs[s]), *prmel_y[s]);
630  }
631  }
632  }
633  }
634  }
635 
636  for (int s=0; s<fes.Size(); ++s)
637  {
638  delete vdofs2[s];
639  delete vdofs[s];
640  delete el_x[s];
641  delete ael_x[s];
642  }
643 
644  for (int s=0; s<paramfes.Size(); ++s)
645  {
646  delete prmvdofs2[s];
647  delete prmvdofs[s];
648  delete prmel_y[s];
649  delete prmel_x[s];
650  }
651 
652 }
653 
654 void ParametricBNLForm::MultBlocked(const BlockVector &bx,
655  const BlockVector &dx,
656  BlockVector &by) const
657 {
658 
659  Array<Array<int> *>vdofs(fes.Size());
660  Array<Array<int> *>vdofs2(fes.Size());
661  Array<Vector *> el_x(fes.Size());
662  Array<const Vector *> el_x_const(fes.Size());
663  Array<Vector *> el_y(fes.Size());
664  Array<const FiniteElement *> fe(fes.Size());
665  Array<const FiniteElement *> fe2(fes.Size());
666 
667  ElementTransformation *T;
668 
669  Array<Array<int> *>prmvdofs(paramfes.Size());
670  Array<Array<int> *>prmvdofs2(paramfes.Size());
671  Array<Vector *> prmel_x(paramfes.Size());
672  Array<const Vector *> prmel_x_const(paramfes.Size());
673  Array<const FiniteElement *> prmfe(paramfes.Size());
674  Array<const FiniteElement *> prmfe2(paramfes.Size());
675 
676  by = 0.0;
677  for (int s=0; s<fes.Size(); ++s)
678  {
679  el_x_const[s] = el_x[s] = new Vector();
680  el_y[s] = new Vector();
681  vdofs[s] = new Array<int>;
682  vdofs2[s] = new Array<int>;
683  }
684 
685  for (int s=0; s<paramfes.Size(); ++s)
686  {
687  prmel_x_const[s] = prmel_x[s] = new Vector();
688  prmvdofs[s] = new Array<int>;
689  prmvdofs2[s] = new Array<int>;
690  }
691 
692  if (dnfi.Size())
693  {
694  for (int i = 0; i < fes[0]->GetNE(); ++i)
695  {
696  T = fes[0]->GetElementTransformation(i);
697  for (int s = 0; s < fes.Size(); ++s)
698  {
699  fes[s]->GetElementVDofs(i, *(vdofs[s]));
700  fe[s] = fes[s]->GetFE(i);
701  bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
702  }
703 
704  for (int s = 0; s < paramfes.Size(); ++s)
705  {
706  paramfes[s]->GetElementVDofs(i, *(prmvdofs[s]));
707  prmfe[s] = paramfes[s]->GetFE(i);
708  dx.GetBlock(s).GetSubVector(*(prmvdofs[s]), *prmel_x[s]);
709  }
710 
711 
712  for (int k = 0; k < dnfi.Size(); ++k)
713  {
714  dnfi[k]->AssembleElementVector(fe,prmfe, *T,
715  el_x_const, prmel_x_const, el_y);
716 
717  for (int s=0; s<fes.Size(); ++s)
718  {
719  if (el_y[s]->Size() == 0) { continue; }
720  by.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
721  }
722  }
723  }
724  }
725 
726  if (fnfi.Size())
727  {
728  Mesh *mesh = fes[0]->GetMesh();
729  FaceElementTransformations *tr;
730 
731  for (int i = 0; i < mesh->GetNumFaces(); ++i)
732  {
733  tr = mesh->GetInteriorFaceTransformations(i);
734  if (tr != NULL)
735  {
736  for (int s=0; s<fes.Size(); ++s)
737  {
738  fe[s] = fes[s]->GetFE(tr->Elem1No);
739  fe2[s] = fes[s]->GetFE(tr->Elem2No);
740 
741  fes[s]->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
742  fes[s]->GetElementVDofs(tr->Elem2No, *(vdofs2[s]));
743 
744  vdofs[s]->Append(*(vdofs2[s]));
745 
746  bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
747  }
748 
749  for (int s=0; s<paramfes.Size(); ++s)
750  {
751  prmfe[s] = paramfes[s]->GetFE(tr->Elem1No);
752  prmfe2[s] = paramfes[s]->GetFE(tr->Elem2No);
753 
754  paramfes[s]->GetElementVDofs(tr->Elem1No, *(prmvdofs[s]));
755  paramfes[s]->GetElementVDofs(tr->Elem2No, *(prmvdofs2[s]));
756 
757  prmvdofs[s]->Append(*(prmvdofs2[s]));
758 
759  dx.GetBlock(s).GetSubVector(*(prmvdofs[s]), *prmel_x[s]);
760  }
761 
762  for (int k = 0; k < fnfi.Size(); ++k)
763  {
764 
765  fnfi[k]->AssembleFaceVector(fe, fe2, prmfe, prmfe2, *tr, el_x_const,
766  prmel_x_const, el_y);
767 
768  for (int s=0; s<fes.Size(); ++s)
769  {
770  if (el_y[s]->Size() == 0) { continue; }
771  by.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
772  }
773  }
774  }
775  }
776  }
777 
778  if (bfnfi.Size())
779  {
780  Mesh *mesh = fes[0]->GetMesh();
781  FaceElementTransformations *tr;
782  // Which boundary attributes need to be processed?
783  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
784  mesh->bdr_attributes.Max() : 0);
785  bdr_attr_marker = 0;
786  for (int k = 0; k < bfnfi.Size(); ++k)
787  {
788  if (bfnfi_marker[k] == NULL)
789  {
790  bdr_attr_marker = 1;
791  break;
792  }
793  Array<int> &bdr_marker = *bfnfi_marker[k];
794  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
795  "invalid boundary marker for boundary face integrator #"
796  << k << ", counting from zero");
797  for (int i = 0; i < bdr_attr_marker.Size(); ++i)
798  {
799  bdr_attr_marker[i] |= bdr_marker[i];
800  }
801  }
802 
803  for (int i = 0; i < mesh->GetNBE(); ++i)
804  {
805  const int bdr_attr = mesh->GetBdrAttribute(i);
806  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
807 
808  tr = mesh->GetBdrFaceTransformations(i);
809  if (tr != NULL)
810  {
811  for (int s=0; s<fes.Size(); ++s)
812  {
813  fe[s] = fes[s]->GetFE(tr->Elem1No);
814  fe2[s] = fes[s]->GetFE(tr->Elem1No);
815 
816  fes[s]->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
817  bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
818  }
819 
820  for (int s=0; s<paramfes.Size(); ++s)
821  {
822  prmfe[s] = paramfes[s]->GetFE(tr->Elem1No);
823  prmfe2[s] = paramfes[s]->GetFE(tr->Elem1No);
824 
825  paramfes[s]->GetElementVDofs(tr->Elem1No, *(prmvdofs[s]));
826  dx.GetBlock(s).GetSubVector(*(prmvdofs[s]), *prmel_x[s]);
827  }
828 
829 
830  for (int k = 0; k < bfnfi.Size(); ++k)
831  {
832  if (bfnfi_marker[k] &&
833  (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
834 
835  bfnfi[k]->AssembleFaceVector(fe, fe2, prmfe, prmfe2, *tr, el_x_const,
836  prmel_x_const, el_y);
837 
838  for (int s=0; s<fes.Size(); ++s)
839  {
840  if (el_y[s]->Size() == 0) { continue; }
841  by.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
842  }
843  }
844  }
845  }
846  }
847 
848  for (int s=0; s<fes.Size(); ++s)
849  {
850  delete vdofs2[s];
851  delete vdofs[s];
852  delete el_y[s];
853  delete el_x[s];
854  }
855 
856  for (int s=0; s<paramfes.Size(); ++s)
857  {
858  delete prmvdofs2[s];
859  delete prmvdofs[s];
860  delete prmel_x[s];
861  }
862 
863 }
864 
865 const BlockVector &ParametricBNLForm::Prolongate(const BlockVector &bx) const
866 {
867  MFEM_VERIFY(bx.Size() == Width(), "invalid input BlockVector size");
868 
869  if (needs_prolongation)
870  {
871  aux1.Update(block_offsets);
872  for (int s = 0; s < fes.Size(); s++)
873  {
874  P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
875  }
876  return aux1;
877  }
878  return bx;
879 }
880 
881 const BlockVector &ParametricBNLForm::ParamProlongate(const BlockVector &bx)
882 const
883 {
884  MFEM_VERIFY(bx.Size() == paramwidth, "invalid input BlockVector size");
885 
886  if (prmneeds_prolongation)
887  {
888  prmaux1.Update(paramblock_offsets);
889  for (int s = 0; s < paramfes.Size(); s++)
890  {
891  Pparam[s]->Mult(bx.GetBlock(s), prmaux1.GetBlock(s));
892  }
893  return prmaux1;
894  }
895  return bx;
896 }
897 
898 
899 void ParametricBNLForm::ParamMult(const Vector &x, Vector &y) const
900 {
901  BlockVector bx(const_cast<Vector&>(x), paramblock_trueOffsets);
902  BlockVector by(y, paramblock_trueOffsets);
903 
904  const BlockVector &pbx = ParamProlongate(bx);
905 
906  if (prmneeds_prolongation)
907  {
908  prmaux2.Update(paramblock_offsets);
909  }
910  BlockVector &pby = prmneeds_prolongation ? prmaux2 : by;
911 
912  xs.Update(const_cast<BlockVector&>(pbx), paramblock_offsets);
913  ys.Update(pby, paramblock_offsets);
914 
915  MultParamBlocked(xsv,adv,xs,ys);
916 
917  for (int s = 0; s < paramfes.Size(); s++)
918  {
919  if (cPparam[s])
920  {
921  cPparam[s]->MultTranspose(pby.GetBlock(s), by.GetBlock(s));
922  }
923  by.GetBlock(s).SetSubVector(*paramess_tdofs[s], 0.0);
924  }
925 }
926 
927 
928 void ParametricBNLForm::Mult(const Vector &x, Vector &y) const
929 {
930 
931  BlockVector bx(const_cast<Vector&>(x), block_trueOffsets);
932  BlockVector by(y, block_trueOffsets);
933 
934  const BlockVector &pbx = Prolongate(bx);
935 
936  if (needs_prolongation)
937  {
938  aux2.Update(block_offsets);
939  }
940  BlockVector &pby = needs_prolongation ? aux2 : by;
941 
942  xs.Update(const_cast<BlockVector&>(pbx), block_offsets);
943  ys.Update(pby, block_offsets);
944  MultBlocked(xs,xdv,ys);
945 
946  for (int s = 0; s < fes.Size(); s++)
947  {
948  if (cP[s])
949  {
950  cP[s]->MultTranspose(pby.GetBlock(s), by.GetBlock(s));
951  }
952  by.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
953  }
954 }
955 
956 void ParametricBNLForm::ComputeGradientBlocked(const BlockVector &bx,
957  const BlockVector &dx) const
958 {
959  const int skip_zeros = 0;
960  Array<Array<int> *> vdofs(fes.Size());
961  Array<Array<int> *> vdofs2(fes.Size());
962  Array<Vector *> el_x(fes.Size());
963  Array<const Vector *> el_x_const(fes.Size());
964  Array2D<DenseMatrix *> elmats(fes.Size(), fes.Size());
965  Array<const FiniteElement *>fe(fes.Size());
966  Array<const FiniteElement *>fe2(fes.Size());
967 
968  ElementTransformation * T;
969 
970  Array<Array<int> *> prmvdofs(paramfes.Size());
971  Array<Array<int> *> prmvdofs2(paramfes.Size());
972  Array<Vector *> prmel_x(paramfes.Size());
973  Array<const Vector *> prmel_x_const(paramfes.Size());
974  Array<const FiniteElement *>prmfe(paramfes.Size());
975  Array<const FiniteElement *>prmfe2(paramfes.Size());
976 
977  for (int i=0; i<fes.Size(); ++i)
978  {
979  el_x_const[i] = el_x[i] = new Vector();
980  vdofs[i] = new Array<int>;
981  vdofs2[i] = new Array<int>;
982  for (int j=0; j<fes.Size(); ++j)
983  {
984  elmats(i,j) = new DenseMatrix();
985  }
986  }
987 
988  for (int i=0; i<fes.Size(); ++i)
989  {
990  for (int j=0; j<fes.Size(); ++j)
991  {
992  if (Grads(i,j) != NULL)
993  {
994  *Grads(i,j) = 0.0;
995  }
996  else
997  {
998  Grads(i,j) = new SparseMatrix(fes[i]->GetVSize(),
999  fes[j]->GetVSize());
1000  }
1001  }
1002  }
1003 
1004  for (int i=0; i<paramfes.Size(); ++i)
1005  {
1006  prmel_x_const[i] = prmel_x[i] = new Vector();
1007  prmvdofs[i] = new Array<int>;
1008  prmvdofs2[i] = new Array<int>;
1009  }
1010 
1011  if (dnfi.Size())
1012  {
1013  for (int i = 0; i < fes[0]->GetNE(); ++i)
1014  {
1015  T = fes[0]->GetElementTransformation(i);
1016 
1017  for (int s = 0; s < fes.Size(); ++s)
1018  {
1019  fe[s] = fes[s]->GetFE(i);
1020  fes[s]->GetElementVDofs(i, *vdofs[s]);
1021  bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
1022  }
1023 
1024  for (int s = 0; s < paramfes.Size(); ++s)
1025  {
1026  prmfe[s] = paramfes[s]->GetFE(i);
1027  paramfes[s]->GetElementVDofs(i, *prmvdofs[s]);
1028  dx.GetBlock(s).GetSubVector(*prmvdofs[s], *prmel_x[s]);
1029  }
1030 
1031  for (int k = 0; k < dnfi.Size(); ++k)
1032  {
1033  dnfi[k]->AssembleElementGrad(fe,prmfe,*T, el_x_const, prmel_x_const, elmats);
1034 
1035  for (int j=0; j<fes.Size(); ++j)
1036  {
1037  for (int l=0; l<fes.Size(); ++l)
1038  {
1039  if (elmats(j,l)->Height() == 0) { continue; }
1040  Grads(j,l)->AddSubMatrix(*vdofs[j], *vdofs[l],
1041  *elmats(j,l), skip_zeros);
1042  }
1043  }
1044  }
1045  }
1046  }
1047 
1048  if (fnfi.Size())
1049  {
1050  FaceElementTransformations *tr;
1051  Mesh *mesh = fes[0]->GetMesh();
1052 
1053  for (int i = 0; i < mesh->GetNumFaces(); ++i)
1054  {
1055  tr = mesh->GetInteriorFaceTransformations(i);
1056 
1057  for (int s=0; s < fes.Size(); ++s)
1058  {
1059  fe[s] = fes[s]->GetFE(tr->Elem1No);
1060  fe2[s] = fes[s]->GetFE(tr->Elem2No);
1061 
1062  fes[s]->GetElementVDofs(tr->Elem1No, *vdofs[s]);
1063  fes[s]->GetElementVDofs(tr->Elem2No, *vdofs2[s]);
1064  vdofs[s]->Append(*(vdofs2[s]));
1065 
1066  bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
1067  }
1068 
1069  for (int s=0; s < paramfes.Size(); ++s)
1070  {
1071  prmfe[s] = paramfes[s]->GetFE(tr->Elem1No);
1072  prmfe2[s] = paramfes[s]->GetFE(tr->Elem2No);
1073 
1074  paramfes[s]->GetElementVDofs(tr->Elem1No, *prmvdofs[s]);
1075  paramfes[s]->GetElementVDofs(tr->Elem2No, *prmvdofs2[s]);
1076  prmvdofs[s]->Append(*(prmvdofs2[s]));
1077 
1078  dx.GetBlock(s).GetSubVector(*prmvdofs[s], *prmel_x[s]);
1079  }
1080 
1081  for (int k = 0; k < fnfi.Size(); ++k)
1082  {
1083  fnfi[k]->AssembleFaceGrad(fe, fe2, prmfe, prmfe2, *tr, el_x_const,
1084  prmel_x_const, elmats);
1085  for (int j=0; j<fes.Size(); ++j)
1086  {
1087  for (int l=0; l<fes.Size(); ++l)
1088  {
1089  if (elmats(j,l)->Height() == 0) { continue; }
1090  Grads(j,l)->AddSubMatrix(*vdofs[j], *vdofs[l],
1091  *elmats(j,l), skip_zeros);
1092  }
1093  }
1094  }
1095  }
1096  }
1097 
1098  if (bfnfi.Size())
1099  {
1100  FaceElementTransformations *tr;
1101  Mesh *mesh = fes[0]->GetMesh();
1102 
1103  // Which boundary attributes need to be processed?
1104  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
1105  mesh->bdr_attributes.Max() : 0);
1106  bdr_attr_marker = 0;
1107  for (int k = 0; k < bfnfi.Size(); ++k)
1108  {
1109  if (bfnfi_marker[k] == NULL)
1110  {
1111  bdr_attr_marker = 1;
1112  break;
1113  }
1114  Array<int> &bdr_marker = *bfnfi_marker[k];
1115  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
1116  "invalid boundary marker for boundary face integrator #"
1117  << k << ", counting from zero");
1118  for (int i = 0; i < bdr_attr_marker.Size(); ++i)
1119  {
1120  bdr_attr_marker[i] |= bdr_marker[i];
1121  }
1122  }
1123 
1124  for (int i = 0; i < mesh->GetNBE(); ++i)
1125  {
1126  const int bdr_attr = mesh->GetBdrAttribute(i);
1127  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
1128 
1129  tr = mesh->GetBdrFaceTransformations(i);
1130  if (tr != NULL)
1131  {
1132  for (int s = 0; s < fes.Size(); ++s)
1133  {
1134  fe[s] = fes[s]->GetFE(tr->Elem1No);
1135  fe2[s] = fe[s];
1136 
1137  fes[s]->GetElementVDofs(i, *vdofs[s]);
1138  bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
1139  }
1140 
1141  for (int s = 0; s < paramfes.Size(); ++s)
1142  {
1143  prmfe[s] = paramfes[s]->GetFE(tr->Elem1No);
1144  prmfe2[s] = prmfe[s];
1145 
1146  paramfes[s]->GetElementVDofs(i, *prmvdofs[s]);
1147  dx.GetBlock(s).GetSubVector(*prmvdofs[s], *prmel_x[s]);
1148  }
1149 
1150 
1151  for (int k = 0; k < bfnfi.Size(); ++k)
1152  {
1153  bfnfi[k]->AssembleFaceGrad(fe, fe2, prmfe, prmfe2, *tr, el_x_const,
1154  prmel_x_const, elmats);
1155  for (int l=0; l<fes.Size(); ++l)
1156  {
1157  for (int j=0; j<fes.Size(); ++j)
1158  {
1159  if (elmats(j,l)->Height() == 0) { continue; }
1160  Grads(j,l)->AddSubMatrix(*vdofs[j], *vdofs[l],
1161  *elmats(j,l), skip_zeros);
1162  }
1163  }
1164  }
1165  }
1166  }
1167  }
1168 
1169  if (!Grads(0,0)->Finalized())
1170  {
1171  for (int i=0; i<fes.Size(); ++i)
1172  {
1173  for (int j=0; j<fes.Size(); ++j)
1174  {
1175  Grads(i,j)->Finalize(skip_zeros);
1176  }
1177  }
1178  }
1179 
1180  for (int i=0; i<fes.Size(); ++i)
1181  {
1182  for (int j=0; j<fes.Size(); ++j)
1183  {
1184  delete elmats(i,j);
1185  }
1186  delete vdofs2[i];
1187  delete vdofs[i];
1188  delete el_x[i];
1189  }
1190 
1191  for (int i=0; i<paramfes.Size(); ++i)
1192  {
1193  delete prmvdofs2[i];
1194  delete prmvdofs[i];
1195  delete prmel_x[i];
1196  }
1197 }
1198 
1199 BlockOperator& ParametricBNLForm::GetGradient(const Vector &x) const
1200 {
1201  BlockVector bx(const_cast<Vector&>(x), block_trueOffsets);
1202  const BlockVector &pbx = Prolongate(bx);
1203 
1204  ComputeGradientBlocked(pbx, xdv);
1205 
1206  Array2D<SparseMatrix *> mGrads(fes.Size(), fes.Size());
1207  mGrads = Grads;
1208  if (needs_prolongation)
1209  {
1210  for (int s1 = 0; s1 < fes.Size(); ++s1)
1211  {
1212  for (int s2 = 0; s2 < fes.Size(); ++s2)
1213  {
1214  delete cGrads(s1, s2);
1215  cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
1216  mGrads(s1, s2) = cGrads(s1, s2);
1217  }
1218  }
1219  }
1220 
1221  for (int s = 0; s < fes.Size(); ++s)
1222  {
1223  for (int i = 0; i < ess_tdofs[s]->Size(); ++i)
1224  {
1225  for (int j = 0; j < fes.Size(); ++j)
1226  {
1227  if (s == j)
1228  {
1229  mGrads(s, s)->EliminateRowCol((*ess_tdofs[s])[i],
1230  Matrix::DIAG_ONE);
1231  }
1232  else
1233  {
1234  mGrads(s, j)->EliminateRow((*ess_tdofs[s])[i]);
1235  mGrads(j, s)->EliminateCol((*ess_tdofs[s])[i]);
1236  }
1237  }
1238  }
1239  }
1240 
1241  delete BlockGrad;
1242  BlockGrad = new BlockOperator(block_trueOffsets);
1243  for (int i = 0; i < fes.Size(); ++i)
1244  {
1245  for (int j = 0; j < fes.Size(); ++j)
1246  {
1247  BlockGrad->SetBlock(i, j, mGrads(i, j));
1248  }
1249  }
1250  return *BlockGrad;
1251 
1252 }
1253 
1254 ParametricBNLForm::~ParametricBNLForm()
1255 {
1256  delete BlockGrad;
1257  for (int i=0; i<fes.Size(); ++i)
1258  {
1259  for (int j=0; j<fes.Size(); ++j)
1260  {
1261  delete Grads(i,j);
1262  delete cGrads(i,j);
1263  }
1264  delete ess_tdofs[i];
1265  }
1266 
1267  for (int i=0; i<paramfes.Size(); ++i)
1268  {
1269  delete paramess_tdofs[i];
1270  }
1271 
1272  for (int i = 0; i < dnfi.Size(); ++i)
1273  {
1274  delete dnfi[i];
1275  }
1276 
1277  for (int i = 0; i < fnfi.Size(); ++i)
1278  {
1279  delete fnfi[i];
1280  }
1281 
1282  for (int i = 0; i < bfnfi.Size(); ++i)
1283  {
1284  delete bfnfi[i];
1285  }
1286 
1287 }
1288 
1289 }
mfem::Array::Size
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
mfem::Vector::SetSubVector
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:551
mfem::Mesh::GetBdrAttribute
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1203
mfem::ParametricBNLForm::SetParamEssentialBC
virtual void SetParamEssentialBC(const Array< Array< int > * > &bdr_attr_is_ess, Array< Vector * > &rhs)
Set the essential boundary conditions on the parametric fields.
Definition: paramnonlinearform.cpp:412
mfem::ParametricBNLForm::P
Array< const Operator * > P
Array of pointers to the prolongation matrix of fes, may be NULL.
Definition: paramnonlinearform.hpp:153
mfem::BlockVector
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
mfem::RAP
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:2485
mfem::ParametricBNLForm::is_serial
bool is_serial
Indicator if the Operator is part of a parallel run.
Definition: paramnonlinearform.hpp:165
mfem::Mesh::GetNBE
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:849
mfem::ParametricBNLForm::fes
Array< FiniteElementSpace * > fes
FE spaces on which the form lives.
Definition: paramnonlinearform.hpp:102
mfem::Operator::Width
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
mfem::Array::Copy
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:851
mfem::FaceElementTransformations
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:467
mfem::Vector::GetSubVector
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:525
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
mfem::BlockVector::Update
void Update(double *data, const Array< int > &bOffsets)
Update method.
Definition: blockvector.cpp:85
mfem::ParametricBNLFormIntegrator::AssembleElementGrad
virtual void AssembleElementGrad(const Array< const FiniteElement * > &el, const Array< const FiniteElement * > &pel, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array< const Vector * > &pelfun, const Array2D< DenseMatrix * > &elmats)
Assemble the local gradient matrix.
Definition: paramnonlinearform.cpp:46
mfem::ParametricBNLForm::SetAdjointFields
virtual void SetAdjointFields(const Vector &av) const
Set the adjoint fields.
Definition: paramnonlinearform.cpp:356
mfem::ParametricBNLFormIntegrator::AssembleFaceGrad
virtual void AssembleFaceGrad(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, const Array< const FiniteElement * > &pel1, const Array< const FiniteElement * > &pel2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array< const Vector * > &pelfun, const Array2D< DenseMatrix * > &elmats)
Assemble the local gradient matrix on faces of the elements.
Definition: paramnonlinearform.cpp:32
mfem::ParametricBNLForm::Pparam
Array< const Operator * > Pparam
Array of pointers to the prolongation matrix of paramfes, may be NULL.
Definition: paramnonlinearform.hpp:156
mfem::ParametricBNLForm::GetEnergy
virtual double GetEnergy(const Vector &x) const
Computes the energy for a state vector x.
Definition: paramnonlinearform.cpp:389
mfem::ParametricBNLFormIntegrator::AssemblePrmFaceVector
virtual void AssemblePrmFaceVector(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, const Array< const FiniteElement * > &pel1, const Array< const FiniteElement * > &pel2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array< const Vector * > &alfun, const Array< const Vector * > &pelfun, const Array< Vector * > &pelvect)
Definition: paramnonlinearform.cpp:98
mfem::Mesh::GetInteriorFaceTransformations
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
Definition: mesh.hpp:1153
mfem::ParametricBNLForm::paramess_tdofs
Array< Array< int > * > paramess_tdofs
Definition: paramnonlinearform.hpp:150
mfem::ParametricBNLForm::paramfes
Array< FiniteElementSpace * > paramfes
FE spaces for the parametric fields.
Definition: paramnonlinearform.hpp:105
mfem::ParametricBNLFormIntegrator::AssemblePrmElementVector
virtual void AssemblePrmElementVector(const Array< const FiniteElement * > &el, const Array< const FiniteElement * > &pel, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array< const Vector * > &alfun, const Array< const Vector * > &pelfun, const Array< Vector * > &pelvec)
Perform the local action on the parameters of the BNLFormIntegrator.
Definition: paramnonlinearform.cpp:85
mfem::ParametricBNLForm::AddBdrFaceIntegrator
void AddBdrFaceIntegrator(ParametricBNLFormIntegrator *nlfi)
Adds new Boundary Face Integrator.
Definition: paramnonlinearform.hpp:252
mfem::Mesh::GetNumFaces
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:4915
mfem::ParametricBNLForm::bfnfi
Array< ParametricBNLFormIntegrator * > bfnfi
Set of Boundary Face Integrators to be assembled (added).
Definition: paramnonlinearform.hpp:117
mfem::ParametricBNLForm::GetEnergyBlocked
double GetEnergyBlocked(const BlockVector &bx, const BlockVector &dx) const
Definition: paramnonlinearform.cpp:257
mfem::SparseMatrix
Data type sparse matrix.
Definition: sparsemat.hpp:41
mfem::Array::Append
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:746
mfem::mfem_error
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:153
mfem::ParametricBNLForm::fnfi
Array< ParametricBNLFormIntegrator * > fnfi
Set of interior face Integrators to be assembled (added).
Definition: paramnonlinearform.hpp:114
mfem::Array::Max
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
mfem::ParametricBNLForm::~ParametricBNLForm
virtual ~ParametricBNLForm()
Destructor.
Definition: paramnonlinearform.cpp:1254
mfem::ParametricBNLForm::MultBlocked
void MultBlocked(const BlockVector &bx, const BlockVector &dx, BlockVector &by) const
Definition: paramnonlinearform.cpp:654
mfem::ParametricBNLForm::Mult
virtual void Mult(const Vector &x, Vector &y) const
Definition: paramnonlinearform.cpp:928
mfem::ParametricBNLForm::prmneeds_prolongation
bool prmneeds_prolongation
Indicator if the Operator needs prolongation on assembly.
Definition: paramnonlinearform.hpp:171
mfem::ParametricBNLForm::SetEssentialBC
virtual void SetEssentialBC(const Array< Array< int > * > &bdr_attr_is_ess, Array< Vector * > &rhs)
Set the essential boundary conditions.
Definition: paramnonlinearform.cpp:395
mfem::Operator::DIAG_ONE
Set the diagonal value to one.
Definition: operator.hpp:49
mfem::Mesh::GetBdrFaceTransformations
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Definition: mesh.cpp:1020
mfem::ParametricBNLForm::SetParamFields
virtual void SetParamFields(const Vector &dv) const
Set the parameters/design fields.
Definition: paramnonlinearform.cpp:373
mfem::ParametricBNLForm::MultParamBlocked
void MultParamBlocked(const BlockVector &bx, const BlockVector &ax, const BlockVector &dx, BlockVector &dy) const
Definition: paramnonlinearform.cpp:429
mfem::Array2D
Dynamic 2D array using row-major layout.
Definition: array.hpp:347
mfem::ParametricBNLForm::cGrads
Array2D< SparseMatrix * > cGrads
Definition: paramnonlinearform.hpp:136
mfem::ParametricBNLForm::cP
Array< const SparseMatrix * > cP
Array of results of dynamic-casting P to SparseMatrix pointer.
Definition: paramnonlinearform.hpp:159
mfem::Mesh::bdr_attributes
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:204
mfem::ParametricBNLForm::SetStateFields
virtual void SetStateFields(const Vector &xv) const
Set the state fields.
Definition: paramnonlinearform.cpp:340
mfem::Array::SetSize
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:674
mfem::Array::PartialSum
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition: array.cpp:103
mfem::ParametricBNLForm::ParamMult
virtual void ParamMult(const Vector &x, Vector &y) const
Definition: paramnonlinearform.cpp:899
mfem::Operator::height
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
mfem::ParametricBNLForm::ess_tdofs
Array< Array< int > * > ess_tdofs
Definition: paramnonlinearform.hpp:147
mfem::ParametricBNLForm::bfnfi_marker
Array< Array< int > * > bfnfi_marker
Definition: paramnonlinearform.hpp:118
mfem::ParametricBNLForm::GetGradient
virtual BlockOperator & GetGradient(const Vector &x) const
Definition: paramnonlinearform.cpp:1199
mfem::ParametricBNLForm::ParametricBNLForm
ParametricBNLForm()
Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
Definition: paramnonlinearform.cpp:114
mfem::ParametricBNLForm::Prolongate
const BlockVector & Prolongate(const BlockVector &bx) const
Definition: paramnonlinearform.cpp:865
mfem::ParametricBNLForm::ComputeGradientBlocked
void ComputeGradientBlocked(const BlockVector &bx, const BlockVector &dx) const
Specialized version of GetGradient() for BlockVector.
Definition: paramnonlinearform.cpp:956
mfem::ParametricBNLForm::Grads
Array2D< SparseMatrix * > Grads
Definition: paramnonlinearform.hpp:136
mfem::ParametricBNLForm::ParamProlongate
const BlockVector & ParamProlongate(const BlockVector &bx) const
Definition: paramnonlinearform.cpp:881
mfem::Vector
Vector data type.
Definition: vector.hpp:60
mfem::ParametricBNLFormIntegrator::AssembleElementVector
virtual void AssembleElementVector(const Array< const FiniteElement * > &el, const Array< const FiniteElement * > &pel, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array< const Vector * > &pelfun, const Array< Vector * > &elvec)
Perform the local action of the BlockNonlinearFormIntegrator.
Definition: paramnonlinearform.cpp:58
mfem::ParametricBNLForm::cPparam
Array< const SparseMatrix * > cPparam
Array of results of dynamic-casting Pparam to SparseMatrix pointer.
Definition: paramnonlinearform.hpp:162
mfem::ParametricBNLFormIntegrator::AssembleFaceVector
virtual void AssembleFaceVector(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, const Array< const FiniteElement * > &pel1, const Array< const FiniteElement * > &pel2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array< const Vector * > &pelfun, const Array< Vector * > &elvect)
Definition: paramnonlinearform.cpp:70
mfem::ParametricBNLForm::dnfi
Array< ParametricBNLFormIntegrator * > dnfi
Set of Domain Integrators to be assembled (added).
Definition: paramnonlinearform.hpp:111
s
RefCoord s[3]
Definition: ncmesh_tables.hpp:158
mfem::BlockOperator
A class to handle Block systems in a matrix-free implementation.
Definition: blockoperator.hpp:34
mfem::ParametricBNLForm::SetSpaces
void SetSpaces(Array< FiniteElementSpace * > &statef, Array< FiniteElementSpace * > &paramf)
(Re)initialize the ParametricBNLForm.
Definition: paramnonlinearform.cpp:130
mfem::Operator::width
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
mfem::BlockOperator::SetBlock
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
Definition: blockoperator.cpp:57
mfem::ParametricBNLForm::needs_prolongation
bool needs_prolongation
Indicator if the Operator needs prolongation on assembly.
Definition: paramnonlinearform.hpp:168
mfem::BlockVector::GetBlock
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:90
mfem::ParametricBNLFormIntegrator::GetElementEnergy
virtual double GetElementEnergy(const Array< const FiniteElement * > &el, const Array< const FiniteElement * > &pel, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array< const Vector * > &pelfun)
Compute the local energy.
Definition: paramnonlinearform.cpp:20
Generated on Thu Jul 29 2021 15:58:41 for MFEM by   doxygen 1.8.5