Geomi
NOXGroup.hpp
1 #ifndef DEF_COMMON_NOXGROUP
2 #define DEF_COMMON_NOXGROUP
3 
4 #include <NOX_Common.H>
5 #include <NOX_Abstract_Group.H>
6 
7 #include <Eigen/Dense>
8 #include <Eigen/Geometry>
9 
10 //#include "Common_NOXVector.hpp"
11 //#include "Common_Abstract_NOXStep.hpp"
12 
13 template <typename T_Q, int T_SIZE_MULTIVECT = 1>
14 class NOXGroup : public virtual NOX::Abstract::Group {
15 public:
16  NOXGroup<T_Q,T_SIZE_MULTIVECT> (Abstract::NOXStep<T_Q,T_SIZE_MULTIVECT>& step)
17  : m_step(step),
18  m_xVector(step.getInitialGuess())
19  { resetIsValid(); }
20 
21  ~NOXGroup<T_Q,T_SIZE_MULTIVECT> ( )
22  { }
23 
24  NOX::Abstract::Group&
25  operator= (const NOXGroup<T_Q,T_SIZE_MULTIVECT>& g)
26  {
27  if (this != &g) {
28  m_xVector = g.m_xVector;
29 
30  m_isValidF = g.m_isValidF;
31  m_isValidJacobian = g.m_isValidJacobian;
32  m_isValidGradient = g.m_isValidGradient;
33  m_isValidNewton = g.m_isValidNewton;
34 
35  if (m_isValidF)
36  m_fVector = g.m_fVector;
37 
38  if (m_isValidJacobian)
39  m_jacobianMatrix = g.m_jacobianMatrix;
40 
41  if (m_isValidGradient)
42  m_gradientVector = g.m_gradientVector;
43 
44  if (m_isValidNewton)
45  m_newtonVector = g.m_newtonVector;
46  }
47 
48  return *this;
49  }
50 
51  NOX::Abstract::Group&
52  operator= (const NOX::Abstract::Group& g)
53  { return operator=(dynamic_cast<const NOXGroup<T_Q,T_SIZE_MULTIVECT>&>(g)); }
54 
55  void
56  setX (const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& x)
57  {
58  resetIsValid();
59  m_xVector = x;
60  }
61 
62  void
63  setX (const NOX::Abstract::Vector& x)
64  { setX(dynamic_cast<const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>&>(x)); }
65 
66  void
67  computeX (const NOXGroup<T_Q,T_SIZE_MULTIVECT>& grp, const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& d, double step)
68  {
69  resetIsValid();
70  m_xVector = grp.m_xVector + step*d;
71  }
72 
73  void
74  computeX (const NOX::Abstract::Group& grp, const NOX::Abstract::Vector& d, double step)
75  {
76  const NOXGroup<T_Q,T_SIZE_MULTIVECT>& g = dynamic_cast<const NOXGroup<T_Q,T_SIZE_MULTIVECT>&>(grp);
77  const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& v = dynamic_cast<const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>&>(d);
78  computeX(g,v,step);
79  }
80 
81  NOX::Abstract::Group::ReturnType
82  computeF ()
83  {
84  if (m_isValidF)
85  return NOX::Abstract::Group::Ok;
86 
87  m_isValidF = m_step.computeF(m_fVector,m_xVector);
88 
89  return (m_isValidF) ? (NOX::Abstract::Group::Ok) : (NOX::Abstract::Group::Failed);
90  }
91 
92  NOX::Abstract::Group::ReturnType
93  computeJacobian ()
94  {
95  if (m_isValidJacobian)
96  return NOX::Abstract::Group::Ok;
97 
98  m_isValidJacobian = m_step.computeJacobian(m_jacobianMatrix,m_xVector);
99 
100  return (m_isValidJacobian) ? (NOX::Abstract::Group::Ok) : (NOX::Abstract::Group::Failed);
101  }
102 
103  NOX::Abstract::Group::ReturnType
104  computeGradient ()
105  {
106  if (m_isValidGradient)
107  return NOX::Abstract::Group::Ok;
108 
109  if (!m_isValidF) {
110  std::cerr
111  << "ERROR: NOXGroup::computeGradient() - F is out of date wrt X"
112  << std::endl;
113  return NOX::Abstract::Group::BadDependency;
114  }
115 
116  if (!m_isValidJacobian) {
117  std::cerr
118  << "ERROR: NOXGroup::computeGradient() - Jacobian is out of date wrt X"
119  << std::endl;
120  return NOX::Abstract::Group::BadDependency;
121  }
122 
123  m_gradientVector = m_jacobianMatrix.transpose()*m_fVector;
124  m_isValidGradient = true;
125 
126  return NOX::Abstract::Group::Ok;
127  }
128 
129  NOX::Abstract::Group::ReturnType
130  computeNewton (Teuchos::ParameterList& p)
131  {
132  if (m_isValidNewton)
133  return NOX::Abstract::Group::Ok;
134 
135  if (!m_isValidF) {
136  std::cerr
137  << "ERROR: NOXGroup::computeNewton() - invalid F"
138  << std::endl;
139  throw "NOX Error";
140  }
141 
142  if (!m_isValidJacobian) {
143  std::cerr
144  << "ERROR: NOXGroup::computeNewton() - invalid Jacobian"
145  << std::endl;
146  throw "NOX Error";
147  }
148 
149  NOX::Abstract::Group::ReturnType status = applyJacobianInverse(p,m_fVector,m_newtonVector);
150  m_isValidNewton = (status == NOX::Abstract::Group::Ok);
151 
152  m_newtonVector = -m_newtonVector;
153 
154  return status;
155  }
156 
157  NOX::Abstract::Group::ReturnType
158  applyJacobian ( const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& input,
159  NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& result) const
160  {
161  if (!m_isValidJacobian)
162  return NOX::Abstract::Group::BadDependency;
163  result = m_jacobianMatrix*input;
164  return NOX::Abstract::Group::Ok;
165  }
166 
167  NOX::Abstract::Group::ReturnType
168  applyJacobian ( const NOX::Abstract::Vector& input,
169  NOX::Abstract::Vector& result) const
170  {
171  const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& v = dynamic_cast<const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>&>(input);
172  NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& res = dynamic_cast<NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>&>(result);
173  return applyJacobian(v,res);
174  }
175 
176  NOX::Abstract::Group::ReturnType
177  applyJacobianTranspose (const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& input,
178  NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& result) const
179  {
180  if (!m_isValidJacobian)
181  return NOX::Abstract::Group::BadDependency;
182  result = m_jacobianMatrix.transpose()*input;
183  return NOX::Abstract::Group::Ok;
184  }
185 
186  NOX::Abstract::Group::ReturnType
187  applyJacobianTranspose (const NOX::Abstract::Vector& input,
188  NOX::Abstract::Vector& result) const
189  {
190  const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& v
191  = dynamic_cast<const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>&>(input);
192  NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& res
193  = dynamic_cast<NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>&>(result);
194  return applyJacobianTranspose(v,res);
195  }
196 
197  NOX::Abstract::Group::ReturnType
198  applyJacobianInverse ( Teuchos::ParameterList& p,
199  const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& input,
200  NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& result) const
201  {
202  if (!m_isValidJacobian) {
203  std::cerr
204  << "ERROR: NOXGroup::applyJacobianInverse() - invalid Jacobian"
205  << std::endl;
206  throw "NOX Error";
207  }
208  // TODO : tester si le jacobien est inversible
209  result = m_jacobianMatrix.inverse()*input;
210  return NOX::Abstract::Group::Ok;
211  }
212 
213  NOX::Abstract::Group::ReturnType
214  applyJacobianInverse ( Teuchos::ParameterList& p,
215  const NOX::Abstract::Vector& input,
216  NOX::Abstract::Vector& result) const
217  {
218  const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& v
219  = dynamic_cast<const NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>&>(input);
220  NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>& res
221  = dynamic_cast<NOXVector<T_Q::DOF*T_SIZE_MULTIVECT>&>(result);
222  return applyJacobianInverse(p,v,res);
223  }
224 
225  // TODO : Overloader applyJacobianInverseMultiVector avec une meilleure implementation que celle par defaut
226  // (qui appelle plusieurs fois applyJacobianInverse)
227 
228  bool
229  isF () const
230  { return m_isValidF; }
231 
232  bool
233  isJacobian () const
234  { return m_isValidJacobian; }
235 
236  bool
237  isGradient () const
238  { return m_isValidGradient; }
239 
240  bool
241  isNewton () const
242  { return m_isValidNewton; }
243 
244  const NOX::Abstract::Vector&
245  getX () const
246  { return m_xVector; }
247 
248  const NOX::Abstract::Vector&
249  getF () const
250  { return m_fVector; }
251 
252  double
253  getNormF () const
254  {
255  if (!m_isValidF) {
256  std::cerr
257  << "ERROR: NOXGroup::getNormF() - invalid F, please call computeF() first"
258  << std::endl;
259  throw "NOX Error";
260  }
261  return m_fVector.norm();
262  }
263 
264  const NOX::Abstract::Vector&
265  getGradient () const
266  { return m_gradientVector; }
267 
268  const NOX::Abstract::Vector&
269  getNewton () const
270  { return m_newtonVector; }
271 
272  Teuchos::RCP<const NOX::Abstract::Vector>
273  getXPtr () const
274  { return Teuchos::RCP<const NOX::Abstract::Vector>(&m_xVector,false); }
275 
276  Teuchos::RCP<const NOX::Abstract::Vector>
277  getFPtr () const
278  { return Teuchos::RCP<const NOX::Abstract::Vector>(&m_fVector,false); }
279 
280  Teuchos::RCP<const NOX::Abstract::Vector>
281  getGradientPtr () const
282  { return Teuchos::RCP<const NOX::Abstract::Vector>(&m_gradientVector,false); }
283 
284  Teuchos::RCP<const NOX::Abstract::Vector>
285  getNewtonPtr () const
286  { return Teuchos::RCP<const NOX::Abstract::Vector>(&m_newtonVector,false); }
287 
288  Teuchos::RCP<NOX::Abstract::Group>
289  clone (NOX::CopyType type = NOX::DeepCopy) const
290  {
291  Teuchos::RCP<NOX::Abstract::Group> tmp;
292  tmp = Teuchos::rcp(new NOXGroup<T_Q,T_SIZE_MULTIVECT>(*this));
293  return tmp;
294  }
295 
296  void
297  print () const
298  {
299  std::cout << "x = " << m_xVector << std::endl;
300  if (m_isValidF) {
301  std::cout << "F(x) = " << m_fVector << std::endl;
302  } else {
303  std::cout << "F(x) has not been computed" << std::endl;
304  }
305  std::cout << std::endl;
306  }
307 
308 protected:
309  void
310  resetIsValid ()
311  {
312  m_isValidF = false;
313  m_isValidJacobian = false;
314  m_isValidGradient = false;
315  m_isValidNewton = false;
316  }
317 
318 protected:
319  NOXVector<T_Q::DOF*T_SIZE_MULTIVECT> m_xVector;
320  NOXVector<T_Q::DOF*T_SIZE_MULTIVECT> m_fVector;
321  NOXVector<T_Q::DOF*T_SIZE_MULTIVECT> m_newtonVector;
322  NOXVector<T_Q::DOF*T_SIZE_MULTIVECT> m_gradientVector;
323  Eigen::Matrix<double,T_Q::DOF*T_SIZE_MULTIVECT,T_Q::DOF*T_SIZE_MULTIVECT> m_jacobianMatrix;
324 
325  Abstract::NOXStep<T_Q,T_SIZE_MULTIVECT>& m_step;
326 
327  bool m_isValidF;
328  bool m_isValidJacobian;
329  bool m_isValidGradient;
330  bool m_isValidNewton;
331 };
332 
333 #endif
NOXGroup
Definition: NOXGroup.hpp:14
Generated by   doxygen 1.8.13