Chombo + EB  3.2
LSProblemImplem.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
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 
34 #include "NamespaceHeader.H"
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_VAL);
98  m_upperBound.assign(HUGE_VAL);
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_VAL;
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 it2 = a_ccm.m_moments.begin();
184  it2 != a_ccm.m_moments.end();++it2)
185  {
186  IvDim mono = it2->first;
187  Real val = it2->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  }
212 
213 
214  m_lowerBound[variableNdx] = lobnd;
215  m_upperBound[variableNdx] = hibnd;
216 
217  }
218 
219 
220 }
221 
222 // Calculate upper and lower bounds for a moment
223 template<int dim> void LSProblem<dim>::momentBounds(Real & a_lobnd,
224  Real & a_hibnd,
225  const IvDim & a_mono,
226  const IFData<dim> & a_IFData)
227 {
228  a_lobnd = 0.0;
229  a_hibnd = 0.0;
230  // Enumerate though every vertex, and integrate the contributions
231  // from the (local) cell center to the vertex in local coordinates
233  // pout() << a_mono << endl;
234  for (cornerIt = a_IFData.m_cornerSigns.begin();
235  cornerIt != a_IFData.m_cornerSigns.end();
236  ++cornerIt)
237  {
238  const typename IFData<dim>::Vertex & vertex = cornerIt->first;
239  // represent the vertex as an RvDim in cell centered coordinates
240  typename IFData<dim>::RvDim corner;
241  for (int idir = 0; idir < dim; ++idir)
242  {
243  corner[idir] = vertex[idir] - 0.5;
244  corner[idir] *= a_IFData.m_cellCenterCoord.m_dx[idir];
245  }
246 
247  // compute coordinates of corner in local coordinates
248  typename IFData<dim>::RvDim cornerCoord
249  = a_IFData.m_localCoord.convert(corner,a_IFData.m_cellCenterCoord);
250 
251  Real partialCellMag = 1.0;
252  for (int idir = 0; idir < dim ; ++idir)
253  {
254  Real pPlus1 =(Real)( a_mono[idir]+1);
255  partialCellMag *= POW(cornerCoord[idir],pPlus1)/pPlus1;
256  if (cornerCoord[idir] < 0.0) partialCellMag *=-1.0;
257  }
258  if (partialCellMag > 0.0) a_hibnd += partialCellMag;
259  if (partialCellMag < 0.0) a_lobnd += partialCellMag;
260  }
261 }
262 
263 template<int dim> LSProblem<dim>::~LSProblem()
264 {
265  if (m_matrix != NULL)
266  {
267  // free m_matrix
268  int numRows = dim*m_numP;
269  int numCols = m_numP + m_numPLess1;
270  freeArray(numRows,numCols,m_matrix);
271 
272  }
273 }
274 
275 // this constructor is used when just a list of monomials is wanted
276 template<int dim> LSProblem<dim>::LSProblem(const int & a_degreeP,
277  const bool & a_useConstraints)
278  :m_degreeP(a_degreeP),
279  m_numActiveBounds(0),
280  m_useConstraints(a_useConstraints)
281 {
286  m_matrix = NULL;
287 }
288 
289 // constructor for solving a LS problem
290 template<int dim> LSProblem<dim>::LSProblem(const int & a_order,
291  const int & a_degreeP,
292  const bool & a_useConstraints,
293  const IndexTM<Real,dim> & a_normal)
294  :m_order(a_order),
295  m_degreeP(a_degreeP),
296  m_numActiveBounds(0),
297  m_useConstraints(a_useConstraints),
298  m_normal(a_normal)
299 {
304  setMatrix();
305 
308 
309  if (a_useConstraints)
310  {
313  }
314 }
315 
316 //
317 template<int dim> int LSProblem<dim>::invertNormalEq(const Vector<Real> & a_rhs,
318  Vector<Real> & a_residual)
319 {
320  m_rhs = a_rhs;
321  int retCode = -1;
322  if (!m_useConstraints)
323  {
324  LSquares lsSolver;
325  lsSolver.LeastSquares(m_matrix,m_unknowns,m_rhs);
326  m_numActiveBounds = 0;
327  retCode = 0;
328  }
329  else
330  {
331  ConstrainedLS cls;
332  ConstrainedLS::LSResult result = cls.solveBoundConstrained(m_unknowns,
333  m_matrix,
334  m_rhs,
335  m_lowerBound,
336  m_upperBound);
337  m_numActiveBounds = cls.numberActiveConstraints();
338  retCode = (result == ConstrainedLS::SUCCESS) ? 0 : -1;
339  switch(result)
340  {
342  break;
344  pout() << "Singular LS problem. Linearly dependent columns (zero normal?)" << endl;
345  break;
347  pout() << "Conflicting upper/lower bounds in LS problem. Zero or bad normal?" << endl;
348  break;
350  MayDay::Error("Malformed LS problem. Underdetermined");
352  pout() << "BVLS failed to converge properly" << endl;
353  break;
354  }
355 
356  }
357 
358  a_residual.resize(3);
359  a_residual[0] = 0.0;
360  a_residual[1] = 0.0;
361  a_residual[2] = 0.0;
362  Real maxRi = 0.0;
363  for (int i = 0 ; i < m_numP*dim ; i++)
364  {
365  Real AtimeX = 0.0;
366  for (int j = 0 ; j < m_numP + m_numPLess1 ; j++)
367  {
368  AtimeX += m_matrix[i][j] * m_unknowns[j];
369  }
370  Real ri = Abs(AtimeX - m_rhs[i]);
371  if (ri > maxRi)
372  {
373  a_residual[0] = ri;
374  maxRi = ri;
375  }
376  a_residual[1] += ri;
377  a_residual[2] += ri * ri;
378  }
379  a_residual[2] = sqrt(a_residual[2]);
380  // pout() << "Residual : " << a_residual[1] << ":" << a_residual[2] << endl;
381  return retCode;
382 }
383 
384 template<int dim> int LSProblem<dim>::factorial(const int & a_n,
385  const int & a_m)
386 {
387  int retval = 1;
388  if (a_n < 0 || a_m < 0)
389  {
390  MayDay::Abort("Attempting n! for n < 0");
391  }
392  for (int i = a_m+1; i <= a_n; ++i)
393  {
394  retval *= i;
395  }
396  return retval;
397 }
398 
399 template<int dim> int LSProblem<dim>::numMonomials(const int & a_monoDegree)
400 {
401  int retval = LARGEINTVAL;
402  if (a_monoDegree == -1)
403  {
404  retval = 0;
405  }
406  else
407  {
408  int bigger;
409  int smaller;
410  if (dim- 1 > a_monoDegree)
411  {
412  bigger = dim-1;
413  smaller = a_monoDegree;
414  }
415  else
416  {
417  smaller = dim-1;
418  bigger = a_monoDegree;
419  }
420  int numerator = factorial(dim - 1 + a_monoDegree,bigger);
421  int denominator = factorial(smaller);
422  // pout() << "numerator = " << numerator << endl;
423  // pout() << "denominator = " << denominator << endl;
424  retval = numerator/denominator;
425  }
426  return retval;
427 }
428 
429 // uses dimension & degree create matrix for overdetermined system
430 template<int dim> void LSProblem<dim>::setMatrix()
431 {
432  // solving m_matrix[x] = b
433 
434  // initializes m_matrix to zeros
435 
436  int numRows = dim*m_numP;
437  int numCols = m_numP + m_numPLess1;
438  allocArray(numRows,numCols,m_matrix);
439 
440  // iterate through the list of mono of Degree P
441  // pout() << " Number of mono of DegreeP = " << m_monoLocP.size() << endl;
442  for (typename PthMomentLoc::const_iterator it = m_monoLocP.begin();
443  it != m_monoLocP.end();++it)
444  {
445  for (int idir = 0; idir< dim; ++idir)
446  {
447  // this entry corresponds to the integral over the boundary
448  int row = dim*(it->second) + idir;
449  int pcol = it->second;
450  m_matrix[row][pcol] = -m_normal[idir];
451 
452  // differentiate the mono and
453  IvDim Dmono;
454  IvDim mono = it->first;
455  int coeff = LARGEINTVAL;
456  // diff(it->first) = coeff*Dmono
457  differentiate(coeff,Dmono,idir,mono);
458  int pLess1Col= LARGEINTVAL;
459  if (coeff == LARGEINTVAL)
460  {
461  MayDay::Abort("problem wth differentiate");
462  }
463  if (coeff > 0)
464  {
465  // find which mono this is in the list
466  if (m_monoLocPLess1.find(Dmono) != m_monoLocPLess1.end())
467  {
468  pLess1Col = m_monoLocPLess1[Dmono] + m_numP;
469  m_matrix[row][pLess1Col] = coeff;
470  }
471  else
472  {
473  MayDay::Abort("can't find derived mono");
474  }
475  }
476  }
477  }
478 }
479 
480 // differentiate a_mono w.r.t. x_idir. Answer = a_coeff*a_Dmono
481 template<int dim> void LSProblem<dim>::differentiate(int & a_coeff,
482  IvDim & a_Dmono,
483  int & a_idir,
484  const IvDim & a_mono)
485 {
486  if (a_mono[a_idir] > 0)
487  {
488  a_coeff = a_mono[a_idir];
489  a_Dmono = a_mono;
490  a_Dmono[a_idir] -= 1;
491  }
492  else if (a_mono[a_idir] == 0)
493  {
494  a_coeff = 0;
495  for (int idir = 0; idir < dim; ++idir)
496  {
497  a_Dmono[idir] = LARGEINTVAL;
498  }
499  }
500  else
501  {
502  MayDay::Abort("Monomial has negative power");
503  }
504 }
505 
506 template<int dim> void LSProblem<dim>::setRhs(const Vector<Real> & a_rhs)
507 {
508 }
509 
510 template<int dim> void LSProblem<dim>::fillMap(PthMomentLoc & a_monoLoc,
511  LocPthMoment & a_locMono,
512  const int & a_degree)
513 {
514  if (a_degree >= 0)
515  {
516  Vector<IvDim> monomials;
517 
518  generateMultiIndices(monomials,a_degree);
519 
520  for (int i = 0; i < monomials.size(); i++)
521  {
522  const IvDim & monomial = monomials[i];
523 
524  a_monoLoc[monomial] = i;
525  a_locMono[i] = monomial;
526  }
527  }
528 }
529 
530 template<int dim> void LSProblem<dim>::allocArray(const int & a_rows,
531  const int & a_cols,
532  Real** & a_A)
533 {
534  a_A = new Real* [a_rows];
535 
536  for (int i = 0; i < a_rows; i++)
537  {
538  a_A[i] = new Real [a_cols];
539 
540  Real* scanA = a_A[i];
541  for (int j = 0; j < a_cols; j++)
542  {
543  *(scanA++) = 0.0;
544  }
545  }
546 }
547 
548 template<int dim> void LSProblem<dim>::freeArray(const int & a_rows,
549  const int & a_cols,
550  Real** & a_A)
551 {
552  for (int i = 0; i < a_rows; i++)
553  {
554  delete[] a_A[i];
555  }
556 
557  delete[] a_A;
558 }
559 
560 template<int dim > void LSProblem<dim>::print(ostream & a_out) const
561 {
562  a_out << "Dim = " << dim << ", degree = " << m_degreeP << '\n';
563  a_out << "m_monoLocP has " << m_monoLocP.size() << " elements" << '\n';
564 
565  for (typename PthMomentLoc::const_iterator it = m_monoLocP.begin();
566  it != m_monoLocP.end();++it)
567  {
568  a_out << "Monomial = " << it->first << ", Loc = " << it->second << '\n';
569  }
570 
571  a_out << "Dim = " << dim << '\n';
572  a_out << "m_locMonoP has " << m_locMonoP.size() << " elements" << '\n';
573 
574  for (typename LocPthMoment::const_iterator it = m_locMonoP.begin();
575  it != m_locMonoP.end();++it)
576  {
577  a_out << "Loc = " << it->first << ", Monomial = " << it->second << '\n';
578  }
579  // one degree lower
580  a_out << "m_locMonoPLess1 has " << m_locMonoPLess1.size() << " elements" << '\n';
581  for (typename PthMomentLoc::const_iterator it = m_monoLocPLess1.begin();
582  it != m_monoLocPLess1.end();++it)
583  {
584  a_out << "Monomial = " << it->first << ", Loc = " << it->second << '\n';
585  }
586 
587  a_out << "m_locMonoPLess1 has " << m_locMonoPLess1.size() << " elements" << '\n';
588 
589  for (typename LocPthMoment::const_iterator it = m_locMonoPLess1.begin();
590  it != m_locMonoPLess1.end();++it)
591  {
592  a_out << "Loc = " << it->first << ", Monomial = " << it->second << '\n';
593  }
594 
595  a_out << "Matrix and rhs for least squares problem of dim = " << dim << endl;
596  outputMatrix();
597  outputRhs();
598  outputUnknowns();
599  outputBounds();
600 }
601 
602 template<int dim> ostream& operator<<(ostream & a_out,
603  LSProblem<dim> & a_lsProblem)
604 {
605  a_lsProblem.print(a_out);
606  return a_out;
607 }
608 
609 template<int dim> void LSProblem<dim>::outputMatrix() const
610 {
611  int rows = m_locMonoP.size()*dim;
612  int cols = m_locMonoPLess1.size() + m_locMonoP.size();
613  pout() << "numRows = " << rows << ", numCols = " << cols << endl;
614  // pout() << "outputting " << name << endl;
615  for (int i = 0; i < rows; i++)
616  {
617  for (int j = 0; j < cols; j++)
618  {
619  pout() << m_matrix[i][j] << " ";
620  }
621  pout() << endl;
622  }
623 }
624 
625 template<int dim> void LSProblem<dim>::outputRhs() const
626 {
627  pout() << "Outputting Rhs" << endl;
628  for (int i = 0; i < m_rhs.size(); i++)
629  {
630  pout() << "m_rhs[" << i << "] = " << m_rhs[i] << endl;
631  }
632 }
633 
634 template<int dim> void LSProblem<dim>::outputUnknowns()const
635 {
636  pout() << "Outputting Unknowns" << endl;
637  for (int i = 0; i < m_unknowns.size(); i++)
638  {
639  pout() << "m_unknowns[" << i << "] = " << m_unknowns[i] << endl;
640  }
641 }
642 
643 template<int dim> void LSProblem<dim>::outputBounds()const
644 {
645  if (m_useConstraints)
646  {
647  pout() << "Outputting Lower/Upper bounds" << endl;
648  for (int i = 0; i < m_unknowns.size(); i++)
649  {
650  pout() << "m_lowerBound[" << setw(2) << i << "] = " << setw(14) << m_lowerBound[i]
651  << " m_upperBound[" << setw(2) << i << "] = " << setw(14) << m_upperBound[i] << endl;
652  }
653  }
654  else
655  {
656  pout() << "Problem is unconstrained" << endl;
657  }
658 }
659 
660 #include "NamespaceFooter.H"
661 
662 #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
void differentiate(int &a_coeff, IvDim &a_Dmono, int &a_idir, const IvDim &a_mono)
Definition: LSProblemImplem.H:481
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:609
int m_numP
Definition: LSProblem.H:186
void outputBounds() const
Definition: LSProblemImplem.H:643
Definition: ConstrainedLS.H:32
int numberActiveConstraints() const
#define LARGEINTVAL
Definition: Notation.H:76
Vector< Real > m_rhs
Definition: LSProblem.H:199
int numMonomials(const int &a_monoDegree)
Definition: LSProblemImplem.H:399
Definition: ConstrainedLS.H:33
Definition: ComputeCutCellMoments.H:35
Real POW(const Real &a_x, const int &a_p)
computes x^p
Definition: Factorial.H:33
void momentBounds(Real &a_lobnd, Real &a_hibnd, const IvDim &a_mono, const IFData< dim > &a_ifData)
Definition: LSProblemImplem.H:223
LSProblem(const int &a_degreeP, const bool &a_useConstraints)
Definition: LSProblemImplem.H:276
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:65
void resize(unsigned int isize)
Definition: Vector.H:346
Definition: ConstrainedLS.H:31
~LSProblem()
Definition: LSProblemImplem.H:263
void fillMap(PthMomentLoc &a_monoLoc, LocPthMoment &a_locMono, const int &a_degree)
Definition: LSProblemImplem.H:510
void print(ostream &a_out) const
Definition: LSProblemImplem.H:560
int factorial(const int &a_n, const int &a_m=0)
Definition: LSProblemImplem.H:384
map< int, IvDim > LocPthMoment
Definition: LSProblem.H:33
Definition: LSquares.H:22
IFData< dim > m_IFData
Definition: CutCellMoments.H:141
Definition: IFData.H:42
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:317
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:506
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.
CoordinateSystem< dim > m_cellCenterCoord
Definition: IFData.H:63
ostream & operator<<(ostream &a_out, LSProblem< dim > &a_lsProblem)
Definition: LSProblemImplem.H:602
int m_numPLess1
Definition: LSProblem.H:193
PthMoment m_moments
Definition: CutCellMoments.H:132
Real ** m_matrix
Definition: LSProblem.H:197
LocPthMoment m_locMonoP
Definition: LSProblem.H:183
T sum() const
Definition: IndexTMI.H:260
CornerSigns m_cornerSigns
Definition: IFData.H:58
map< IvDim, int > PthMomentLoc
Definition: LSProblem.H:32
Real factorial(const int n)
Calculates factorial for an integer.
Definition: Factorial.H:25
size_t size() const
Definition: Vector.H:192
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:135
void allocArray(const int &a_rows, const int &a_cols, Real **&a_A)
Definition: LSProblemImplem.H:530
Definition: ConstrainedLS.H:19
void outputRhs() const
Definition: LSProblemImplem.H:625
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:548
int dim
Definition: EBInterface.H:146
void setMatrix()
Definition: LSProblemImplem.H:430
Definition: CutCellMoments.H:32
Definition: ConstrainedLS.H:34
Vector< Real > m_unknowns
Definition: LSProblem.H:198
void outputUnknowns() const
Definition: LSProblemImplem.H:634
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).