Chombo + EB  3.0
LSProblemImplem.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
8  */
9 #endif
10
11 #ifndef _LSPROBLEMIMPLEM_H_
12 #define _LSPROBLEMIMPLEM_H_
13
14 #if defined(CH_Darwin) && defined(__GNUC__) && ( __GNUC__ == 3 )
15 // deal with the broken isnan()/isinf() in GCC on MacOS
16 #include <unistd.h>
17 #define _GLIBCPP_USE_C99 1
18 #endif
19
20 #include <cmath>
21 #include <cstdio>
22 #include <cstdlib>
23 #include <fstream>
24 #include <iostream>
25 #include <string>
26 using std::endl;
27
28 #include "ConstrainedLS.H"
29 #include "LSquares.H"
30 #include "MultiIndex.H"
31 #include "IFData.H"
32 #include "CutCellMoments.H"
33
35
36 template<int dim> int LSProblem<dim>::nChooseR(int a_n,
37  int a_r)
38 {
39  if (a_r == 0) return 1;
40  int num = 1;
41  int den = 1;
42  for (int i = 1; i <= a_r; ++i)
43  {
44  num *= (a_n+1-i);
45  den *= i;
46  }
47  return num/den;
48 }
49
50 template<int dim> void LSProblem<dim>::monoMaxMin(Real & a_maxVal,
51  Real & a_minVal,
52  const IndexTM<int,dim> & a_mono,
53  const IFData<dim> & a_IFData)
54 {
55  a_maxVal = 0.0;
56  a_minVal = 0.0;
57
58  for (typename IFData<dim>::CornerSigns::const_iterator cornerIt = a_IFData.m_cornerSigns.begin();
59  cornerIt != a_IFData.m_cornerSigns.end();
60  ++cornerIt)
61  {
62  const typename IFData<dim>::Vertex & vertex = cornerIt->first;
63
64  // represent the vertex as an RvDim in cell centered coordinates
65  typename IFData<dim>::RvDim corner;
66  for (int idir = 0; idir < dim; ++idir)
67  {
68  corner[idir] = vertex[idir] - 0.5;
69  corner[idir] *= a_IFData.m_cellCenterCoord.m_dx[idir];
70  }
71
72  // compute coordinates of corner in local coordinates
73  typename IFData<dim>::RvDim cornerCoord = a_IFData.m_localCoord.convert(corner,a_IFData.m_cellCenterCoord);
74
75  // monomial value at the corner
76  Real monoCorner = 1.0;
77  for (int idir = 0; idir < dim ; ++idir)
78  {
79  monoCorner *= pow(cornerCoord[idir],a_mono[idir]);
80  }
81
82  if (monoCorner > a_maxVal)
83  {
84  a_maxVal = monoCorner;
85  }
86
87  if (monoCorner < a_minVal)
88  {
89  a_minVal = monoCorner;
90  }
91  }
92 }
93
94 template<int dim> void LSProblem<dim>::computeBounds(const IndexTM<Real,dim> & a_dx,
95  const CutCellMoments<dim> & a_ccm)
96 {
97  m_lowerBound.assign(-HUGE);
98  m_upperBound.assign(HUGE);
99
100  int nEBBounds = 0;
101  if (m_degreeP == 0 && dim == 2) // || m_degreeP == 1)
102  {
103  nEBBounds = m_numP; // A lower bound for each EB. Haven't figured out a useful u.b. yet
104  }
105
106  // int nVar=m_numP + m_numPLess1; // total variables in problem
107
108  // This constraint is a lower bound on the zero monomial of the EB.
109  // It is based on inverting integral(mono)<Max(mono)*integral(1)
110  if (m_degreeP == 0 && dim == 2)
111  {
112  Real lobnd = 0.0; // minimum, but we can do better
113  IvDim whichMono; // for tracking
114  for (typename CutCellMoments<dim>::PthMoment::const_iterator it = a_ccm.m_EBmoments.begin();
115  it != a_ccm.m_EBmoments.end();++it)
116  {
117  IvDim mono = it->first;
118  Real maxMono = 0.0;
119  Real minMono = 0.0;
120  monoMaxMin(maxMono,minMono,mono,a_ccm.m_IFData);
121  Real val = it->second;
122  if (minMono < 0.0)
123  {
124  Real newLoBnd = val/minMono;
125  if (newLoBnd > lobnd)
126  {
127  lobnd = newLoBnd;
128  whichMono = mono;
129  }
130  }
131  if (maxMono > 0.0)
132  {
133  Real newLoBnd = val/maxMono;
134  if (newLoBnd > lobnd)
135  {
136  lobnd = newLoBnd;
137  whichMono = mono;
138  }
139  }
140  }
141
142  // The constraints matrix fits CiT x + ci >=0
143  // so the rows are variables and the columns are constraints (may be a FORTRAN holdover)
144  IvDim mono = IvDim::Zero;
145
146  int ndx = 0; // todo: guarantee there should be only one for (0,0)
147  int variableNdx = ndx;
148  m_lowerBound[variableNdx]=lobnd;
149  m_upperBound[variableNdx]=HUGE;
150
151  }
152
153  // Now create constraints for the volume integrals.
154  for (typename PthMomentLoc::const_iterator it = m_monoLocPLess1.begin();
155  it != m_monoLocPLess1.end();++it)
156  {
157  IvDim mono = it->first;
158  // The constraints matrix fits CiT x + ci >=0
159  // so the rows are variables and the columns are constraints (may be a FORTRAN holdover)
160  int ndx = it->second;
161  int variableNdx = ndx + m_numP;
162  Real lobnd, hibnd;
163
164  // Compute the lower and upper bounds
165  // based on quadrature of all the regions that contribute
166  // negatively or positively.
167  momentBounds(lobnd,hibnd,mono, a_ccm.m_IFData);
168
169
170
171  bool isZeroDegree = (mono.sum() == 0);
172
173  if (isZeroDegree )
174  {
175  Real zeroDegreeLoBnd = 0.0;
176
177  // Refine the lower bound of the zero degree monomial
178  // based by inverting the identity
179  // integral(mono)/integral(zeromono) < Max(mono)
180  // Note that we already have integral(mono) from previous steps
181  IvDim whichMono;
182  // Real valAtWhichMono=LARGEREALVAL;
183  for (typename CutCellMoments<dim>::PthMoment::const_iterator it = a_ccm.m_moments.begin();
184  it != a_ccm.m_moments.end();++it)
185  {
186  IvDim mono = it->first;
187  Real val = it->second;
188  Real maxMono = 0.0;
189  Real minMono = 0.0;
190  monoMaxMin(maxMono,minMono,mono,a_ccm.m_IFData);
191  if (minMono < 0.0)
192  {
193  Real newLoBnd = val/minMono;
194  if (newLoBnd > zeroDegreeLoBnd)
195  {
196  zeroDegreeLoBnd = newLoBnd;
197  whichMono = mono;
198  }
199  }
200  if (maxMono > 0.0)
201  {
202  Real newLoBnd = val/maxMono;
203  if (newLoBnd > zeroDegreeLoBnd)
204  {
205  zeroDegreeLoBnd = newLoBnd;
206  whichMono = mono;
207  }
208  }
209  }
210  lobnd = Max(lobnd, zeroDegreeLoBnd);
211 #ifdef CH_REALM
212  // This represents the assumption that kappa3d < kappa2d
214
215  //Refine lower bound based on assumption that kappa3d < kappa2d*dz
217  it != a_ccm.m_bdCutCellMoments.end();
218  ++it)
219  {
220  Iv2 bdId = it->first;
221  int idir = bdId[BDID_DIR];
222  EBorVol doVol = VolMoment;
223  Real volSide = it->second.getVol(doVol);
224  if (volSide > largestVol[idir])
225  {
226  largestVol[idir] = volSide;
227  }
228  }
229  Real noUndercutBnd = LARGEREALVAL;
230  if (a_ccm.m_IFData.m_globalCoord.m_dx[0] != a_ccm.m_IFData.m_globalCoord.m_dx[dim-1])
231  {
232  Real vol= largestVol[dim-1];
233  Real normalDx = a_ccm.m_IFData.m_globalCoord.m_dx[dim-1];
234  noUndercutBnd = Min(noUndercutBnd, vol*normalDx);
235  }
236  hibnd = Min(hibnd, noUndercutBnd);
237 #endif
238  }
239
240
241  m_lowerBound[variableNdx] = lobnd;
242  m_upperBound[variableNdx] = hibnd;
243
244  }
245
246
247 }
248
249 // Calculate upper and lower bounds for a moment
250 template<int dim> void LSProblem<dim>::momentBounds(Real & a_lobnd,
251  Real & a_hibnd,
252  const IvDim & a_mono,
253  const IFData<dim> & a_IFData)
254 {
255  a_lobnd = 0.0;
256  a_hibnd = 0.0;
257  // Enumerate though every vertex, and integrate the contributions
258  // from the (local) cell center to the vertex in local coordinates
260  // pout() << a_mono << endl;
261  for (cornerIt = a_IFData.m_cornerSigns.begin();
262  cornerIt != a_IFData.m_cornerSigns.end();
263  ++cornerIt)
264  {
265  const typename IFData<dim>::Vertex & vertex = cornerIt->first;
266  // represent the vertex as an RvDim in cell centered coordinates
267  typename IFData<dim>::RvDim corner;
268  for (int idir = 0; idir < dim; ++idir)
269  {
270  corner[idir] = vertex[idir] - 0.5;
271  corner[idir] *= a_IFData.m_cellCenterCoord.m_dx[idir];
272  }
273
274  // compute coordinates of corner in local coordinates
275  typename IFData<dim>::RvDim cornerCoord
276  = a_IFData.m_localCoord.convert(corner,a_IFData.m_cellCenterCoord);
277
278  Real partialCellMag = 1.0;
279  for (int idir = 0; idir < dim ; ++idir)
280  {
281  Real pPlus1 =(Real)( a_mono[idir]+1);
282  partialCellMag *= pow(cornerCoord[idir],pPlus1)/pPlus1;
283  if (cornerCoord[idir] < 0.0) partialCellMag *=-1.0;
284  }
285  if (partialCellMag > 0.0) a_hibnd += partialCellMag;
286  if (partialCellMag < 0.0) a_lobnd += partialCellMag;
287  }
288 }
289
290 template<int dim> LSProblem<dim>::~LSProblem()
291 {
292  if (m_matrix != NULL)
293  {
294  // free m_matrix
295  int numRows = dim*m_numP;
296  int numCols = m_numP + m_numPLess1;
297  freeArray(numRows,numCols,m_matrix);
298
299  }
300 }
301
302 // this constructor is used when just a list of monomials is wanted
303 template<int dim> LSProblem<dim>::LSProblem(const int & a_degreeP,
304  const bool & a_useConstraints)
305  :m_degreeP(a_degreeP),
306  m_numActiveBounds(0),
307  m_useConstraints(a_useConstraints)
308 {
313  m_matrix = NULL;
314 }
315
316 // constructor for solving a LS problem
317 template<int dim> LSProblem<dim>::LSProblem(const int & a_order,
318  const int & a_degreeP,
319  const bool & a_useConstraints,
320  const IndexTM<Real,dim> & a_normal)
321  :m_order(a_order),
322  m_degreeP(a_degreeP),
323  m_numActiveBounds(0),
324  m_useConstraints(a_useConstraints),
325  m_normal(a_normal)
326 {
331  setMatrix();
332
335
336  if (a_useConstraints)
337  {
340  }
341 }
342
343 //
344 template<int dim> int LSProblem<dim>::invertNormalEq(const Vector<Real> & a_rhs,
345  Vector<Real> & a_residual)
346 {
347  m_rhs = a_rhs;
348  int retCode = -1;
349  if (!m_useConstraints)
350  {
351  LSquares lsSolver;
352  lsSolver.LeastSquares(m_matrix,m_unknowns,m_rhs);
353  m_numActiveBounds = 0;
354  retCode = 0;
355  }
356  else
357  {
358  ConstrainedLS cls;
359  ConstrainedLS::LSResult result = cls.solveBoundConstrained(m_unknowns,
360  m_matrix,
361  m_rhs,
362  m_lowerBound,
363  m_upperBound);
364  m_numActiveBounds = cls.numberActiveConstraints();
365  retCode = (result == ConstrainedLS::SUCCESS) ? 0 : -1;
366  switch(result)
367  {
369  break;
371  pout() << "Singular LS problem. Linearly dependent columns (zero normal?)" << endl;
372  break;
374  pout() << "Conflicting upper/lower bounds in LS problem. Zero or bad normal?" << endl;
375  break;
377  MayDay::Error("Malformed LS problem. Underdetermined");
379  pout() << "BVLS failed to converge properly" << endl;
380  break;
381  }
382
383  }
384
385  a_residual.resize(3);
386  a_residual[0] = 0.0;
387  a_residual[1] = 0.0;
388  a_residual[2] = 0.0;
389  Real maxRi = 0.0;
390  for (int i = 0 ; i < m_numP*dim ; i++)
391  {
392  Real AtimeX = 0.0;
393  for (int j = 0 ; j < m_numP + m_numPLess1 ; j++)
394  {
395  AtimeX += m_matrix[i][j] * m_unknowns[j];
396  }
397  Real ri = Abs(AtimeX - m_rhs[i]);
398  if (ri > maxRi)
399  {
400  a_residual[0] = ri;
401  maxRi = ri;
402  }
403  a_residual[1] += ri;
404  a_residual[2] += ri * ri;
405  }
406  a_residual[2] = sqrt(a_residual[2]);
407  // pout() << "Residual : " << a_residual[1] << ":" << a_residual[2] << endl;
408  return retCode;
409 }
410
411 template<int dim> int LSProblem<dim>::factorial(const int & a_n,
412  const int & a_m)
413 {
414  int retval = 1;
415  if (a_n < 0 || a_m < 0)
416  {
417  MayDay::Abort("Attempting n! for n < 0");
418  }
419  for (int i = a_m+1; i <= a_n; ++i)
420  {
421  retval *= i;
422  }
423  return retval;
424 }
425
426 template<int dim> int LSProblem<dim>::numMonomials(const int & a_monoDegree)
427 {
428  int retval = LARGEINTVAL;
429  if (a_monoDegree == -1)
430  {
431  retval = 0;
432  }
433  else
434  {
435  int bigger;
436  int smaller;
437  if (dim- 1 > a_monoDegree)
438  {
439  bigger = dim-1;
440  smaller = a_monoDegree;
441  }
442  else
443  {
444  smaller = dim-1;
445  bigger = a_monoDegree;
446  }
447  int numerator = factorial(dim - 1 + a_monoDegree,bigger);
448  int denominator = factorial(smaller);
449  // pout() << "numerator = " << numerator << endl;
450  // pout() << "denominator = " << denominator << endl;
451  retval = numerator/denominator;
452  }
453  return retval;
454 }
455
456 // uses dimension & degree create matrix for overdetermined system
457 template<int dim> void LSProblem<dim>::setMatrix()
458 {
459  // solving m_matrix[x] = b
460
461  // initializes m_matrix to zeros
462
463  int numRows = dim*m_numP;
464  int numCols = m_numP + m_numPLess1;
465  allocArray(numRows,numCols,m_matrix);
466
467  // iterate through the list of mono of Degree P
468  // pout() << " Number of mono of DegreeP = " << m_monoLocP.size() << endl;
469  for (typename PthMomentLoc::const_iterator it = m_monoLocP.begin();
470  it != m_monoLocP.end();++it)
471  {
472  for (int idir = 0; idir< dim; ++idir)
473  {
474  // this entry corresponds to the integral over the boundary
475  int row = dim*(it->second) + idir;
476  int pcol = it->second;
477  m_matrix[row][pcol] = -m_normal[idir];
478
479  // differentiate the mono and
480  IvDim Dmono;
481  IvDim mono = it->first;
482  int coeff = LARGEINTVAL;
483  // diff(it->first) = coeff*Dmono
484  differentiate(coeff,Dmono,idir,mono);
485  int pLess1Col= LARGEINTVAL;
486  if (coeff == LARGEINTVAL)
487  {
488  MayDay::Abort("problem wth differentiate");
489  }
490  if (coeff > 0)
491  {
492  // find which mono this is in the list
493  if (m_monoLocPLess1.find(Dmono) != m_monoLocPLess1.end())
494  {
495  pLess1Col = m_monoLocPLess1[Dmono] + m_numP;
496  m_matrix[row][pLess1Col] = coeff;
497  }
498  else
499  {
500  MayDay::Abort("can't find derived mono");
501  }
502  }
503  }
504  }
505 }
506
507 // differentiate a_mono w.r.t. x_idir. Answer = a_coeff*a_Dmono
508 template<int dim> void LSProblem<dim>::differentiate(int & a_coeff,
509  IvDim & a_Dmono,
510  int & a_idir,
511  const IvDim & a_mono)
512 {
513  if (a_mono[a_idir] > 0)
514  {
515  a_coeff = a_mono[a_idir];
516  a_Dmono = a_mono;
517  a_Dmono[a_idir] -= 1;
518  }
519  else if (a_mono[a_idir] == 0)
520  {
521  a_coeff = 0;
522  for (int idir = 0; idir < dim; ++idir)
523  {
524  a_Dmono[idir] = LARGEINTVAL;
525  }
526  }
527  else
528  {
529  MayDay::Abort("Monomial has negative power");
530  }
531 }
532
533 template<int dim> void LSProblem<dim>::setRhs(const Vector<Real> & a_rhs)
534 {
535 }
536
537 template<int dim> void LSProblem<dim>::fillMap(PthMomentLoc & a_monoLoc,
538  LocPthMoment & a_locMono,
539  const int & a_degree)
540 {
541  if (a_degree >= 0)
542  {
543  Vector<IvDim> monomials;
544
545  generateMultiIndices(monomials,a_degree);
546
547  for (int i = 0; i < monomials.size(); i++)
548  {
549  const IvDim & monomial = monomials[i];
550
551  a_monoLoc[monomial] = i;
552  a_locMono[i] = monomial;
553  }
554  }
555 }
556
557 template<int dim> void LSProblem<dim>::allocArray(const int & a_rows,
558  const int & a_cols,
559  Real** & a_A)
560 {
561  a_A = new Real* [a_rows];
562
563  for (int i = 0; i < a_rows; i++)
564  {
565  a_A[i] = new Real [a_cols];
566
567  Real* scanA = a_A[i];
568  for (int j = 0; j < a_cols; j++)
569  {
570  *(scanA++) = 0.0;
571  }
572  }
573 }
574
575 template<int dim> void LSProblem<dim>::freeArray(const int & a_rows,
576  const int & a_cols,
577  Real** & a_A)
578 {
579  for (int i = 0; i < a_rows; i++)
580  {
581  delete[] a_A[i];
582  }
583
584  delete[] a_A;
585 }
586
587 template<int dim > void LSProblem<dim>::print(ostream & a_out) const
588 {
589  a_out << "Dim = " << dim << ", degree = " << m_degreeP << '\n';
590  a_out << "m_monoLocP has " << m_monoLocP.size() << " elements" << '\n';
591
592  for (typename PthMomentLoc::const_iterator it = m_monoLocP.begin();
593  it != m_monoLocP.end();++it)
594  {
595  a_out << "Monomial = " << it->first << ", Loc = " << it->second << '\n';
596  }
597
598  a_out << "Dim = " << dim << '\n';
599  a_out << "m_locMonoP has " << m_locMonoP.size() << " elements" << '\n';
600
601  for (typename LocPthMoment::const_iterator it = m_locMonoP.begin();
602  it != m_locMonoP.end();++it)
603  {
604  a_out << "Loc = " << it->first << ", Monomial = " << it->second << '\n';
605  }
606  // one degree lower
607  a_out << "m_locMonoPLess1 has " << m_locMonoPLess1.size() << " elements" << '\n';
608  for (typename PthMomentLoc::const_iterator it = m_monoLocPLess1.begin();
609  it != m_monoLocPLess1.end();++it)
610  {
611  a_out << "Monomial = " << it->first << ", Loc = " << it->second << '\n';
612  }
613
614  a_out << "m_locMonoPLess1 has " << m_locMonoPLess1.size() << " elements" << '\n';
615
616  for (typename LocPthMoment::const_iterator it = m_locMonoPLess1.begin();
617  it != m_locMonoPLess1.end();++it)
618  {
619  a_out << "Loc = " << it->first << ", Monomial = " << it->second << '\n';
620  }
621
622  a_out << "Matrix and rhs for least squares problem of dim = " << dim << endl;
623  outputMatrix();
624  outputRhs();
625  outputUnknowns();
626  outputBounds();
627 }
628
629 template<int dim> ostream& operator<<(ostream & a_out,
630  LSProblem<dim> & a_lsProblem)
631 {
632  a_lsProblem.print(a_out);
633  return a_out;
634 }
635
636 template<int dim> void LSProblem<dim>::outputMatrix() const
637 {
638  int rows = m_locMonoP.size()*dim;
639  int cols = m_locMonoPLess1.size() + m_locMonoP.size();
640  pout() << "numRows = " << rows << ", numCols = " << cols << endl;
641  // pout() << "outputting " << name << endl;
642  for (int i = 0; i < rows; i++)
643  {
644  for (int j = 0; j < cols; j++)
645  {
646  pout() << m_matrix[i][j] << " ";
647  }
648  pout() << endl;
649  }
650 }
651
652 template<int dim> void LSProblem<dim>::outputRhs() const
653 {
654  pout() << "Outputting Rhs" << endl;
655  for (int i = 0; i < m_rhs.size(); i++)
656  {
657  pout() << "m_rhs[" << i << "] = " << m_rhs[i] << endl;
658  }
659 }
660
661 template<int dim> void LSProblem<dim>::outputUnknowns()const
662 {
663  pout() << "Outputting Unknowns" << endl;
664  for (int i = 0; i < m_unknowns.size(); i++)
665  {
666  pout() << "m_unknowns[" << i << "] = " << m_unknowns[i] << endl;
667  }
668 }
669
670 template<int dim> void LSProblem<dim>::outputBounds()const
671 {
672  if (m_useConstraints)
673  {
674  pout() << "Outputting Lower/Upper bounds" << endl;
675  for (int i = 0; i < m_unknowns.size(); i++)
676  {
677  pout() << "m_lowerBound[" << setw(2) << i << "] = " << setw(14) << m_lowerBound[i]
678  << " m_upperBound[" << setw(2) << i << "] = " << setw(14) << m_upperBound[i] << endl;
679  }
680  }
681  else
682  {
683  pout() << "Problem is unconstrained" << endl;
684  }
685 }
686
687 #include "NamespaceFooter.H"
688
689 #endif
std::ostream & pout()
Use this in place of std::cout for program output.
void computeBounds(const IndexTM< Real, dim > &a_dx, const CutCellMoments< dim > &a_ccm)
Definition: LSProblemImplem.H:94
BdCutCellMoments m_bdCutCellMoments
Definition: CutCellMoments.H:144
void differentiate(int &a_coeff, IvDim &a_Dmono, int &a_idir, const IvDim &a_mono)
Definition: LSProblemImplem.H:508
map< IvDim, int, LexLT< IvDim > > PthMomentLoc
Definition: LSProblem.H:32
LSResult solveBoundConstrained(Vector< Real > &a_x, Real **a_A, const Vector< Real > &a_rhs, const Vector< Real > &a_lowerBound, const Vector< Real > &a_upperBound)
void outputMatrix() const
Definition: LSProblemImplem.H:636
int m_numP
Definition: LSProblem.H:186
void outputBounds() const
Definition: LSProblemImplem.H:670
Definition: ConstrainedLS.H:32
int numberActiveConstraints() const
#define LARGEREALVAL
Definition: Notation.H:77
#define LARGEINTVAL
Definition: Notation.H:76
EBorVol
Definition: Notation.H:101
Vector< Real > m_rhs
Definition: LSProblem.H:199
int numMonomials(const int &a_monoDegree)
Definition: LSProblemImplem.H:426
Definition: ConstrainedLS.H:33
Definition: ComputeCutCellMoments.H:35
void momentBounds(Real &a_lobnd, Real &a_hibnd, const IvDim &a_mono, const IFData< dim > &a_ifData)
Definition: LSProblemImplem.H:250
LSProblem(const int &a_degreeP, const bool &a_useConstraints)
Definition: LSProblemImplem.H:303
Definition: ConstrainedLS.H:35
PthMomentLoc m_monoLocPLess1
Definition: LSProblem.H:189
LSResult
Definition: ConstrainedLS.H:29
Definition: IndexTM.H:36
CoordinateSystem< dim > m_localCoord
Definition: IFData.H:58
void resize(unsigned int isize)
Definition: Vector.H:323
Definition: ConstrainedLS.H:31
~LSProblem()
Definition: LSProblemImplem.H:290
void fillMap(PthMomentLoc &a_monoLoc, LocPthMoment &a_locMono, const int &a_degree)
Definition: LSProblemImplem.H:537
void print(ostream &a_out) const
Definition: LSProblemImplem.H:587
int factorial(const int &a_n, const int &a_m=0)
Definition: LSProblemImplem.H:411
map< int, IvDim > LocPthMoment
Definition: LSProblem.H:33
Definition: LSquares.H:22
IFData< dim > m_IFData
Definition: CutCellMoments.H:147
Definition: IFData.H:35
#define BDID_DIR
Definition: Notation.H:88
int nChooseR(int a_n, int a_r)
Definition: LSProblemImplem.H:36
double Real
Definition: REAL.H:33
int invertNormalEq(const Vector< Real > &a_rhs, Vector< Real > &a_residual)
Definition: LSProblemImplem.H:344
T Abs(const T &a_a)
Definition: Misc.H:53
void monoMaxMin(Real &a_maxVal, Real &a_minVal, const IndexTM< int, dim > &a_mono, const IFData< dim > &a_IFData)
Definition: LSProblemImplem.H:50
PthMomentLoc m_monoLocP
Definition: LSProblem.H:182
void setRhs(const Vector< Real > &a_rhs)
Definition: LSProblemImplem.H:533
static void Error(const char *const a_msg=m_nullString, int m_exitCode=CH_DEFAULT_ERROR_CODE)
Print out message to cerr and exit with the specified exit code.
Definition: Notation.H:104
CoordinateSystem< dim > m_cellCenterCoord
Definition: IFData.H:56
ostream & operator<<(ostream &a_out, LSProblem< dim > &a_lsProblem)
Definition: LSProblemImplem.H:629
int m_numPLess1
Definition: LSProblem.H:193
PthMoment m_moments
Definition: CutCellMoments.H:138
Real ** m_matrix
Definition: LSProblem.H:197
LocPthMoment m_locMonoP
Definition: LSProblem.H:183
T sum() const
Definition: IndexTMI.H:260
T Min(const T &a_a, const T &a_b)
Definition: Misc.H:26
CornerSigns m_cornerSigns
Definition: IFData.H:51
size_t size() const
Definition: Vector.H:177
int m_degreeP
Definition: LSProblem.H:170
Vector< Real > m_upperBound
Definition: LSProblem.H:202
T Max(const T &a_a, const T &a_b)
Definition: Misc.H:39
PthMoment m_EBmoments
Definition: CutCellMoments.H:141
void allocArray(const int &a_rows, const int &a_cols, Real **&a_A)
Definition: LSProblemImplem.H:557
Definition: ConstrainedLS.H:19
void outputRhs() const
Definition: LSProblemImplem.H:652
void generateMultiIndices(Vector< IndexTM< int, dim > > &a_indices, const int &a_magnitude)
Definition: MultiIndexImplem.H:18
void LeastSquares(Real **A, Vector< Real > &x, const Vector< Real > &rhs)
void freeArray(const int &a_rows, const int &a_cols, Real **&a_A)
Definition: LSProblemImplem.H:575
int dim
Definition: EBInterface.H:146
void setMatrix()
Definition: LSProblemImplem.H:457
Definition: CutCellMoments.H:32
Definition: ConstrainedLS.H:34
Vector< Real > m_unknowns
Definition: LSProblem.H:198
void outputUnknowns() const
Definition: LSProblemImplem.H:661
Vector< Real > m_lowerBound
Definition: LSProblem.H:201
LocPthMoment m_locMonoPLess1
Definition: LSProblem.H:190
static void Abort(const char *const a_msg=m_nullString)
Print out message to cerr and exit via abort() (if serial) or MPI_Abort() (if parallel).