rsvs3D  0.0.0
Codes for the c++ implementation of the 3D RSVS
RSVScalc_core.cpp
1 #include <Eigen>
2 #include <cmath>
3 #include <cstdlib>
4 #include <fstream>
5 #include <iostream>
6 
7 #include "RSVScalc.hpp"
8 #include "RSVScalctools.hpp"
9 #include "RSVSmath.hpp"
10 #include "matrixtools.hpp"
11 #include "rsvsutils.hpp"
12 #include "snake.hpp"
13 #include "triangulate.hpp"
14 #include "vectorarray.hpp"
15 #include "warning.hpp"
16 
17 using namespace std;
18 using namespace Eigen;
19 using namespace rsvs3d::constants;
20 
21 /*
22 Implementation of the core calculation functions of
23 RSVScalc, this is done to reduce the size of the objects and
24 keep compilation time manageable.
25 
26 */
27 
28 //==========================================
29 // Core class functions
30 //==========================================
31 
32 void RSVScalc::CalcTriangle(const triangle &triIn, const triangulation &triRSVS, bool isObj, bool isConstr,
33  bool isDeriv)
34 {
35  int ii, jj, nj, nCellTarg;
36  int nDvAct;
37  SurfCentroid centreCalc;
38  Volume VolumeCalc;
39  Area AreaCalc;
40  vector<int> subTempVec;
41  HashedVector<int, int> dvListMap;
42  grid::coordlist veccoord;
43  MatrixXd HPos, dPos;
44  MatrixXd HVal, dVal;
45  MatrixXd dConstrPart, HConstrPart, HObjPart;
46  MatrixXd dObjPart;
47  double constrPart, objPart;
48  ArrayVec<double> *HValpnt = NULL, *dValpnt = NULL;
49  double *retVal;
50 
51  auto clock1 = rsvs3d::TimeStamp(NULL, 0);
52  // Prepare calculation steps
53  veccoord = TrianglePointerCoordinates(triIn, triRSVS);
54  dvListMap = TriangleActiveDesignVariables(triIn, triRSVS, this->dvMap, this->useSurfCentreDeriv);
55  TrianglePositionDerivatives(triIn, triRSVS, dvListMap, dPos, HPos, this->useSurfCentreDeriv,
56  this->useSurfCentreHessian);
57 
58  auto clock2 = rsvs3d::TimeStamp(NULL, clock1);
59  this->timer1 += (clock2 - clock1);
60  clock1 = clock2;
61  // Total
62 
63  // Constr and objective
64  nDvAct = dvListMap.vec.size();
65  VolumeCalc.assign(veccoord[0], veccoord[1], veccoord[2]);
66  AreaCalc.assign(veccoord[0], veccoord[1], veccoord[2]);
67 
68  HVal.setZero(9, 9);
69  dVal.setZero(1, 9);
70  dConstrPart.setZero(1, nDvAct);
71  HConstrPart.setZero(nDvAct, nDvAct);
72 
73  if (isConstr)
74  {
75  VolumeCalc.Calc();
76  VolumeCalc.ReturnDatPoint(&retVal, &dValpnt, &HValpnt);
77  ArrayVec2MatrixXd(*HValpnt, HVal);
78  ArrayVec2MatrixXd(*dValpnt, dVal);
79  Deriv1stChainScalar(dVal, dPos, dConstrPart);
80  Deriv2ndChainScalar(dVal, dPos, HVal, HPos, HConstrPart);
81  // HConstrPart = dConstrPart.transpose()*dConstrPart;
82 
83  constrPart = *retVal;
84  // Assign Constraint
85  // and constraint derivative
86  // and Hessian
87  nCellTarg = triIn.connec.celltarg.size();
88  for (ii = 0; ii < nCellTarg; ++ii)
89  {
90  subTempVec = this->constrMap.findall(triIn.connec.celltarg[ii]);
91  nj = subTempVec.size();
92  for (jj = 0; jj < nj; ++jj)
93  {
94  if (subTempVec[jj] != __notfound)
95  {
96  this->constr[subTempVec[jj]] += triIn.connec.constrinfluence[ii] * constrPart;
97  if (isDeriv)
98  {
99  if (this->UseFullMath())
100  {
101  AssignConstraintDerivativesFullMath(*this, triIn, dvListMap, dConstrPart, HConstrPart,
102  subTempVec[jj], ii);
103  }
104  else
105  {
106  AssignConstraintDerivativesSparseMath(*this, triIn, dvListMap, dConstrPart, HConstrPart,
107  subTempVec[jj], ii);
108  }
109  }
110  }
111  else
112  {
113  this->falseaccess++;
114  }
115  }
116  }
117  }
118  clock2 = rsvs3d::TimeStamp(NULL, clock1);
119  this->timer2 += (clock2 - clock1);
120  clock1 = clock2;
121 
122  HVal.setZero(9, 9);
123  dVal.setZero(1, 9);
124  dObjPart.setZero(1, nDvAct);
125  HObjPart.setZero(nDvAct, nDvAct);
126 
127  if (isObj)
128  {
129  AreaCalc.Calc();
130  AreaCalc.ReturnDatPoint(&retVal, &dValpnt, &HValpnt);
131  ArrayVec2MatrixXd(*HValpnt, HVal);
132  ArrayVec2MatrixXd(*dValpnt, dVal);
133  Deriv1stChainScalar(dVal, dPos, dObjPart);
134  Deriv2ndChainScalar(dVal, dPos, HVal, HPos, HObjPart);
135  // HObjPart = dObjPart.transpose()*dObjPart;
136  objPart = *retVal;
137 
138  // Assign to main part of the object
139  // assign objective function
140  this->obj += objPart;
141  // Assign objective derivative
142  if (isDeriv)
143  {
144  for (ii = 0; ii < nDvAct; ++ii)
145  {
146  this->dObj[this->dvMap.find(dvListMap.vec[ii])] += dObjPart(0, ii);
147  if (this->UseFullMath())
148  {
149  for (jj = 0; jj < nDvAct; ++jj)
150  {
151  this->HObj(dvMap.find(dvListMap.vec[jj]), this->dvMap.find(dvListMap.vec[ii])) +=
152  HObjPart(jj, ii);
153  }
154  }
155  else
156  {
157  for (jj = 0; jj < nDvAct; ++jj)
158  {
159  if (!rsvs3d::utils::IsAproxEqual(HObjPart(jj, ii), 0.0))
160  {
161  this->HObj_sparse.coeffRef(dvMap.find(dvListMap.vec[jj]),
162  this->dvMap.find(dvListMap.vec[ii])) += HObjPart(jj, ii);
163  }
164  }
165  }
166  }
167  }
168  }
169  clock2 = rsvs3d::TimeStamp(NULL, clock1);
170  this->timer3 += (clock2 - clock1);
171  clock1 = clock2;
172  // Update active lists of design variables
173  for (ii = 0; ii < nDvAct; ++ii)
174  {
175  this->isDvAct.at(dvMap.find(dvListMap.vec[ii])) = true;
176  }
177  if (nDvAct > 0)
178  {
179  nCellTarg = triIn.connec.celltarg.size();
180  for (ii = 0; ii < nCellTarg; ++ii)
181  {
182  subTempVec = constrMap.findall(triIn.connec.celltarg[ii]);
183  nj = subTempVec.size();
184  for (jj = 0; jj < nj; ++jj)
185  {
186  if (subTempVec[jj] != __notfound)
187  {
188  this->isConstrAct.at(subTempVec[jj]) = true;
189  }
190  }
191  }
192  }
193 
194  // Assign Objective Hessian
195 
196  // Assign Constraint Hessian
197 }
198 
199 void RSVScalc::CalcTriangleFD(const triangle &triIn, const triangulation &triRSVS, bool isObj, bool isConstr,
200  bool isDeriv)
201 {
202  /*Same as calctriangle but the volume derivative is calculated using a
203  Finite Difference.
204 
205  helps isolate errors in RSVSmath.
206  */
207 
208  int ii, jj, nj, nCellTarg;
209  int nDvAct;
210  SurfCentroid centreCalc;
211  Volume VolumeCalc;
212  Area AreaCalc;
213  vector<int> dvList, subTempVec;
214  HashedVector<int, int> dvListMap;
215  grid::coordlist veccoord;
216  MatrixXd HPos, dPos;
217  MatrixXd HVal, dVal;
218  MatrixXd dConstrPart, HConstrPart, HObjPart;
219  MatrixXd dObjPart;
220  double constrPart, objPart;
221  ArrayVec<double> *HValpnt = NULL, *dValpnt = NULL;
222  double *retVal;
223 
224  // Prepare calculation steps
225  veccoord = TrianglePointerCoordinates(triIn, triRSVS);
226  dvListMap = TriangleActiveDesignVariables(triIn, triRSVS, this->dvMap, this->useSurfCentreDeriv);
227  TrianglePositionDerivatives(triIn, triRSVS, dvListMap, dPos, HPos, this->useSurfCentreDeriv,
228  this->useSurfCentreHessian);
229 
230  // Constr and objective
231  nDvAct = dvListMap.vec.size();
232  VolumeCalc.assign(veccoord[0], veccoord[1], veccoord[2]);
233  AreaCalc.assign(veccoord[0], veccoord[1], veccoord[2]);
234 
235  HVal.setZero(9, 9);
236  dVal.setZero(1, 9);
237  dConstrPart.setZero(1, nDvAct);
238  HConstrPart.setZero(nDvAct, nDvAct);
239 
240  if (isConstr)
241  {
242  VolumeCalc.CalcFD();
243  VolumeCalc.ReturnDatPoint(&retVal, &dValpnt, &HValpnt);
244  ArrayVec2MatrixXd(*HValpnt, HVal);
245  ArrayVec2MatrixXd(*dValpnt, dVal);
246  Deriv1stChainScalar(dVal, dPos, dConstrPart);
247  Deriv2ndChainScalar(dVal, dPos, HVal, HPos, HConstrPart);
248  // HConstrPart = dConstrPart.transpose()*dConstrPart;
249  constrPart = *retVal;
250  // Assign Constraint
251  // and constraint derivative
252  // and Hessian
253  nCellTarg = triIn.connec.celltarg.size();
254  for (ii = 0; ii < nCellTarg; ++ii)
255  {
256  subTempVec = this->constrMap.findall(triIn.connec.celltarg[ii]);
257  nj = subTempVec.size();
258  for (jj = 0; jj < nj; ++jj)
259  {
260  if (subTempVec[jj] != __notfound)
261  {
262  this->constr[subTempVec[jj]] += triIn.connec.constrinfluence[ii] * constrPart;
263  if (isDeriv)
264  {
265  if (this->UseFullMath())
266  {
267  AssignConstraintDerivativesFullMath(*this, triIn, dvListMap, dConstrPart, HConstrPart,
268  subTempVec[jj], ii);
269  }
270  else
271  {
272  AssignConstraintDerivativesSparseMath(*this, triIn, dvListMap, dConstrPart, HConstrPart,
273  subTempVec[jj], ii);
274  }
275  }
276  }
277  else
278  {
279  this->falseaccess++;
280  }
281  }
282  }
283  }
284 
285  HVal.setZero(9, 9);
286  dVal.setZero(1, 9);
287  dObjPart.setZero(1, nDvAct);
288  HObjPart.setZero(nDvAct, nDvAct);
289 
290  if (isObj)
291  {
292  AreaCalc.Calc();
293  AreaCalc.ReturnDatPoint(&retVal, &dValpnt, &HValpnt);
294  ArrayVec2MatrixXd(*HValpnt, HVal);
295  ArrayVec2MatrixXd(*dValpnt, dVal);
296  Deriv1stChainScalar(dVal, dPos, dObjPart);
297  Deriv2ndChainScalar(dVal, dPos, HVal, HPos, HObjPart);
298  // HObjPart = dObjPart.transpose()*dObjPart;
299  objPart = *retVal;
300 
301  // Assign to main part of the object
302  // assign objective function
303  this->obj += objPart;
304  // Assign objective derivative
305  if (isDeriv)
306  {
307  for (ii = 0; ii < nDvAct; ++ii)
308  {
309  this->dObj[this->dvMap.find(dvListMap.vec[ii])] += dObjPart(0, ii);
310  if (this->UseFullMath())
311  {
312  for (jj = 0; jj < nDvAct; ++jj)
313  {
314  this->HObj(dvMap.find(dvListMap.vec[jj]), this->dvMap.find(dvListMap.vec[ii])) +=
315  HObjPart(jj, ii);
316  }
317  }
318  else
319  {
320  for (jj = 0; jj < nDvAct; ++jj)
321  {
322  this->HObj_sparse.coeffRef(dvMap.find(dvListMap.vec[jj]),
323  this->dvMap.find(dvListMap.vec[ii])) += HObjPart(jj, ii);
324  }
325  }
326  }
327  }
328  }
329  // Update active lists of design variables
330  for (ii = 0; ii < nDvAct; ++ii)
331  {
332  this->isDvAct.at(dvMap.find(dvListMap.vec[ii])) = true;
333  }
334  if (nDvAct > 0)
335  {
336  nCellTarg = triIn.connec.celltarg.size();
337  for (ii = 0; ii < nCellTarg; ++ii)
338  {
339  subTempVec = constrMap.findall(triIn.connec.celltarg[ii]);
340  nj = subTempVec.size();
341  for (jj = 0; jj < nj; ++jj)
342  {
343  if (subTempVec[jj] != __notfound)
344  {
345  this->isConstrAct.at(subTempVec[jj]) = true;
346  }
347  }
348  }
349  }
350  // Assign Objective Hessian
351 
352  // Assign Constraint Hessian
353 }
354 
355 void RSVScalc::CalcTriangleDirectVolume(const triangle &triIn, const triangulation &triRSVS, bool isObj, bool isConstr,
356  bool isDeriv)
357 {
358  /*Same as calctriangle but the volume derivative is calculated without
359  intermediate steps
360 
361  helps isolate errors on dPos.
362  */
363 
364  int ii, ni, jj, nj, nCellTarg;
365  int subTemp, nDvAct;
366  SurfCentroid centreCalc;
367  Volume2 VolumeCalc2;
368  Volume VolumeCalc;
369  Area AreaCalc;
370  vector<int> dvList, subTempVec, dvOrder;
371  HashedVector<int, int> dvListMap;
372  grid::coordlist veccoord, veccoordvol;
373  vector<double> dvec;
374  MatrixXd HPos, dPos;
375  MatrixXd HVal, dVal;
376  MatrixXd dConstrPart, HConstrPart, HObjPart;
377  MatrixXd dObjPart;
378  double constrPart, objPart;
379  ArrayVec<double> *HValpnt = NULL, *dValpnt = NULL;
380  ArrayVec<double> *HValpnt2 = NULL, *dValpnt2 = NULL;
381  double *retVal = NULL;
382  double *retVal2 = NULL;
383 
384  // Prepare calculation steps
385  veccoord = TrianglePointerCoordinates(triIn, triRSVS);
386  dvListMap = TriangleActiveDesignVariables(triIn, triRSVS, this->dvMap, this->useSurfCentreDeriv);
387  TrianglePositionDerivatives(triIn, triRSVS, dvListMap, dPos, HPos, this->useSurfCentreDeriv,
388  this->useSurfCentreHessian);
389 
390  // Constr and objective
391  nDvAct = dvListMap.vec.size();
392  VolumeCalc.assign(veccoord[0], veccoord[1], veccoord[2]);
393  AreaCalc.assign(veccoord[0], veccoord[1], veccoord[2]);
394 
395  HVal.setZero(9, 9);
396  dVal.setZero(1, 9);
397  dConstrPart.setZero(1, nDvAct);
398  HConstrPart.setZero(nDvAct, nDvAct);
399 
400  // Volume calc assignement
401  dvOrder.assign(3, 0);
402  ni = 3;
403  for (ii = 0; ii < ni; ++ii)
404  {
405  if (triIn.pointtype[ii] == meshtypes::mesh)
406  {
407  dvec[ii] = 0;
408  veccoordvol[1 + ii] = &(triRSVS.meshDep->verts.isearch(triIn.pointind[ii])->coord);
409  veccoordvol[1 + 3 + ii] = &(triRSVS.meshDep->verts.isearch(triIn.pointind[ii])->coord);
410  }
411  else if (triIn.pointtype[ii] == meshtypes::snake)
412  {
413  subTemp = triRSVS.snakeDep->snakeconn.verts.find(triIn.pointind[ii]);
414  dvOrder[ii] = triIn.pointind[ii];
415  dvec[ii] = triRSVS.snakeDep->snaxs(subTemp)->d;
416  veccoordvol[1 + ii] = &(triRSVS.meshDep->verts.isearch(triRSVS.snakeDep->snaxs(subTemp)->fromvert)->coord);
417  veccoordvol[1 + ii + 3] =
418  &(triRSVS.meshDep->verts.isearch(triRSVS.snakeDep->snaxs(subTemp)->tovert)->coord);
419  }
420  else if (triIn.pointtype[ii] == meshtypes::triangulation)
421  {
422  dvec[ii] = 0; // dvec used as a pointer
423  veccoordvol[1 + ii] = (triRSVS.trivert.isearch(triIn.pointind[ii])->coord.retPtr());
424  veccoordvol[1 + ii + 3] = (triRSVS.trivert.isearch(triIn.pointind[ii])->coord.retPtr());
425  }
426  }
427 
428  // Constr and objective
429  VolumeCalc2.assign(veccoordvol);
430  VolumeCalc.assign(veccoord[0], veccoord[1],
431  veccoord[2]);
432  AreaCalc.assign(veccoord[0], veccoord[1], veccoord[2]);
433 
434  HVal.setZero(9, 9);
435  dVal.setZero(1, 9);
436  dConstrPart.setZero(1, nDvAct);
437  HConstrPart.setZero(nDvAct, nDvAct);
438 
439  if (isConstr)
440  {
441  VolumeCalc2.Calc();
442  VolumeCalc.Calc();
443  VolumeCalc.ReturnDatPoint(&retVal, &dValpnt, &HValpnt);
444  VolumeCalc2.ReturnDatPoint(&retVal2, &dValpnt2, &HValpnt2);
445  ArrayVec2MatrixXd(*HValpnt2, HConstrPart);
446  ArrayVec2MatrixXd(*dValpnt2, dConstrPart);
447 
448  if (fabs(*retVal - *retVal2) > 1e-10)
449  {
450  this->falseaccess++;
451  }
452  constrPart = *retVal2;
453 
454  if (dvListMap.find(1044) != __notfound && isDeriv)
455  {
456  cout << endl;
457  DisplayVector(dvOrder);
458  cout << endl;
459  PrintMatrix(dConstrPart);
460  cout << endl;
461  }
462 
463  // Assign Constraint
464  // and constraint derivative
465  // and Hessian
466  nCellTarg = triIn.connec.celltarg.size();
467  for (ii = 0; ii < nCellTarg; ++ii)
468  {
469  subTempVec = this->constrMap.findall(triIn.connec.celltarg[ii]);
470  nj = subTempVec.size();
471  for (jj = 0; jj < nj; ++jj)
472  {
473  if (subTempVec[jj] != __notfound)
474  {
475  this->constr[subTempVec[jj]] += triIn.connec.constrinfluence[ii] * constrPart;
476  if (isDeriv)
477  {
478  if (this->UseFullMath())
479  {
480  AssignConstraintDerivativesFullMath(*this, triIn, dvListMap, dConstrPart, HConstrPart,
481  subTempVec[jj], ii);
482  }
483  else
484  {
485  AssignConstraintDerivativesSparseMath(*this, triIn, dvListMap, dConstrPart, HConstrPart,
486  subTempVec[jj], ii);
487  }
488  }
489  }
490  else
491  {
492  this->falseaccess++;
493  }
494  }
495  }
496  }
497 
498  HVal.setZero(9, 9);
499  dVal.setZero(1, 9);
500  dObjPart.setZero(1, nDvAct);
501  HObjPart.setZero(nDvAct, nDvAct);
502 
503  if (isObj)
504  {
505  AreaCalc.Calc();
506  AreaCalc.ReturnDatPoint(&retVal, &dValpnt, &HValpnt);
507  ArrayVec2MatrixXd(*HValpnt, HVal);
508  ArrayVec2MatrixXd(*dValpnt, dVal);
509  Deriv1stChainScalar(dVal, dPos, dObjPart);
510  Deriv2ndChainScalar(dVal, dPos, HVal, HPos, HObjPart);
511  objPart = *retVal;
512 
513  // Assign to main part of the object
514  // assign objective function
515  this->obj += objPart;
516  // Assign objective derivative
517  if (isDeriv)
518  {
519  for (ii = 0; ii < nDvAct; ++ii)
520  {
521  this->dObj[this->dvMap.find(dvListMap.vec[ii])] += dObjPart(0, ii);
522  if (this->UseFullMath())
523  {
524  for (jj = 0; jj < nDvAct; ++jj)
525  {
526  this->HObj(dvMap.find(dvListMap.vec[jj]), this->dvMap.find(dvListMap.vec[ii])) +=
527  HObjPart(jj, ii);
528  }
529  }
530  else
531  {
532  for (jj = 0; jj < nDvAct; ++jj)
533  {
534  this->HObj_sparse.coeffRef(dvMap.find(dvListMap.vec[jj]),
535  this->dvMap.find(dvListMap.vec[ii])) += HObjPart(jj, ii);
536  }
537  }
538  }
539  }
540  }
541  // Update active lists of design variables
542  for (ii = 0; ii < nDvAct; ++ii)
543  {
544  this->isDvAct.at(dvMap.find(dvListMap.vec[ii])) = true;
545  }
546  if (nDvAct > 0)
547  {
548  nCellTarg = triIn.connec.celltarg.size();
549  for (ii = 0; ii < nCellTarg; ++ii)
550  {
551  subTempVec = constrMap.findall(triIn.connec.celltarg[ii]);
552  nj = subTempVec.size();
553  for (jj = 0; jj < nj; ++jj)
554  {
555  if (subTempVec[jj] != __notfound)
556  {
557  this->isConstrAct.at(subTempVec[jj]) = true;
558  }
559  }
560  }
561  }
562 
563  // Assign Objective Hessian
564 
565  // Assign Constraint Hessian
566 }
567 
568 // =====================================
569 //
570 //
571 // Flag change is temporary
572 #pragma GCC diagnostic push
573 #pragma GCC diagnostic ignored "-Wunused-parameter"
574 void RSVScalc::CalcTriangleEdgeLength(const triangle &triIn, const triangulation &triRSVS, bool isObj, bool isConstr,
575  bool isDeriv)
576 {
577 #pragma GCC diagnostic pop
578 
579  int ii, ni, jj, nj, nCellTarg;
580  int subTemp, subTemp1, subTemp2, subTemp3, nDvAct;
581  SurfCentroid centreCalc;
582  LengthEdge EdgeCalc;
583  Area AreaCalc;
584  vector<int> dvList, subTempVec;
585  HashedVector<int, int> dvListMap;
586  grid::coordlist veccoord;
587  MatrixXd HPos, dPos;
588  MatrixXd HVal, dVal, HVal2, dVal2;
589  MatrixXd dConstrPart, HConstrPart, HObjPart;
590  MatrixXd dObjPart;
591  double objPart;
592  ArrayVec<double> *HValpnt = NULL, *dValpnt = NULL;
593  double *retVal;
594 
595  veccoord.reserve(3);
596  ni = 3;
597  for (ii = 0; ii < ni; ++ii)
598  {
599  if (triIn.pointtype[ii] == meshtypes::mesh)
600  {
601  veccoord.push_back(&(triRSVS.meshDep->verts.isearch(triIn.pointind[ii])->coord));
602  }
603  else if (triIn.pointtype[ii] == meshtypes::snake)
604  {
605  veccoord.push_back(&(triRSVS.snakeDep->snakeconn.verts.isearch(triIn.pointind[ii])->coord));
606  }
607  else if (triIn.pointtype[ii] == meshtypes::triangulation)
608  {
609  veccoord.push_back((triRSVS.trivert.isearch(triIn.pointind[ii])->coord.retPtr()));
610  }
611  }
612 
613  // Constr and objective
614  // CoordFunc.assign(int pRepI,vector<double> &pRep);
615  AreaCalc.assign(veccoord[0], veccoord[1], veccoord[2]);
616 
617  // Active DV lists
618 
619  for (ii = 0; ii < ni; ++ii)
620  {
621  if (triIn.pointtype[ii] == meshtypes::snake && this->dvMap.find(triIn.pointind[ii]) != __notfound)
622  {
623  dvList.push_back(triIn.pointind[ii]);
624  }
625  else if (triIn.pointtype[ii] == meshtypes::triangulation && false)
626  {
627  subTemp = triRSVS.trivert.find(triIn.pointind[ii]);
628  nj = triRSVS.trisurf(subTemp)->indvert.size();
629  for (jj = 0; jj < nj; ++jj)
630  {
631  if (triRSVS.trisurf(subTemp)->typevert[jj] == meshtypes::snake)
632  {
633  dvList.push_back(triRSVS.trisurf(subTemp)->indvert[jj]);
634  }
635  }
636  }
637  }
638  dvListMap.vec = dvList;
639  sort(dvListMap.vec);
640  unique(dvListMap.vec);
641  dvListMap.GenerateHash();
642  nDvAct = dvListMap.vec.size();
643 
644  // Positional Derivatives
645 
646  // HERE -> function to calculate SurfCentroid (dc/dd)^T Hm (dc/dd)
647 
648  HPos.setZero(nDvAct, nDvAct);
649  dPos.setZero(9, nDvAct);
650  for (ii = 0; ii < ni; ++ii)
651  {
652  if (triIn.pointtype[ii] == meshtypes::snake && this->dvMap.find(triIn.pointind[ii]) != __notfound)
653  {
654  subTemp = triRSVS.snakeDep->snaxs.find(triIn.pointind[ii]);
655  subTemp1 = triRSVS.meshDep->verts.find(triRSVS.snakeDep->snaxs(subTemp)->fromvert);
656  subTemp2 = triRSVS.meshDep->verts.find(triRSVS.snakeDep->snaxs(subTemp)->tovert);
657  subTemp3 = dvListMap.find(triIn.pointind[ii]);
658  for (jj = 0; jj < 3; ++jj)
659  {
660  dPos(ii * 3 + jj, subTemp3) +=
661  triRSVS.meshDep->verts(subTemp2)->coord[jj] - triRSVS.meshDep->verts(subTemp1)->coord[jj];
662  }
663  }
664  }
665  // Total
666 
667  if (false)
668  {
669  RSVS3D_ERROR_ARGUMENT("2D Area constraint not implemented yet");
670  HVal.setZero(9, 9);
671  dVal.setZero(1, 9);
672  dConstrPart.setZero(1, nDvAct);
673  HConstrPart.setZero(nDvAct, nDvAct);
674  AreaCalc.Calc();
675  AreaCalc.ReturnDatPoint(&retVal, &dValpnt, &HValpnt);
676  ArrayVec2MatrixXd(*HValpnt, HVal);
677  ArrayVec2MatrixXd(*dValpnt, dVal);
678  Deriv1stChainScalar(dVal, dPos, dConstrPart);
679  Deriv2ndChainScalar(dVal, dPos, HVal, HPos, HConstrPart);
680  }
681 
682  if (isObj)
683  {
684  HVal.setZero(9, 9);
685  dVal.setZero(1, 9);
686  objPart = 0;
687  for (int ed = 0; ed < 3; ++ed)
688  {
689  if (triIn.pointtype[ed] == meshtypes::snake && triIn.pointtype[(ed + 1) % 3] == meshtypes::snake)
690  {
691  HVal2.setZero(6, 6);
692  dVal2.setZero(1, 6);
693  EdgeCalc.assign(0, (triRSVS.snakeDep->snakeconn.verts.isearch(triIn.pointind[ed])->coord));
694  EdgeCalc.assign(1, (triRSVS.snakeDep->snakeconn.verts.isearch(triIn.pointind[(ed + 1) % 3])->coord));
695  EdgeCalc.Calc();
696  EdgeCalc.ReturnDatPoint(&retVal, &dValpnt, &HValpnt);
697  objPart += *retVal;
698  ArrayVec2MatrixXd(*HValpnt, HVal2);
699  ArrayVec2MatrixXd(*dValpnt, dVal2);
700  for (ii = 0; ii < 6; ++ii)
701  {
702  dVal(0, (ed * 3 + ii) % 9) += dVal2(0, ii);
703  for (jj = 0; jj < 6; ++jj)
704  {
705  HVal((ed * 3 + jj) % 9, (ed * 3 + ii) % 9) += HVal2(0, jj);
706  }
707  }
708  }
709  }
710 
711  Deriv1stChainScalar(dVal, dPos, dObjPart);
712  Deriv2ndChainScalar(dVal, dPos, HVal, HPos, HObjPart);
713  // Assign to main part of the object
714  // assign objective function
715  this->obj += objPart;
716  // Assign objective derivative
717  if (isDeriv)
718  {
719  for (ii = 0; ii < nDvAct; ++ii)
720  {
721  this->dObj[this->dvMap.find(dvListMap.vec[ii])] += dObjPart(0, ii);
722  for (jj = 0; jj < nDvAct; ++jj)
723  {
724  this->HObj(dvMap.find(dvListMap.vec[jj]), this->dvMap.find(dvListMap.vec[ii])) += HObjPart(jj, ii);
725  }
726  }
727  }
728  }
729 
730  // Update active lists of design variables
731  for (ii = 0; ii < nDvAct; ++ii)
732  {
733  this->isDvAct.at(dvMap.find(dvListMap.vec[ii])) = true;
734  }
735  nCellTarg = triIn.connec.celltarg.size();
736  for (ii = 0; ii < nCellTarg; ++ii)
737  {
738  subTempVec = constrMap.findall(triIn.connec.celltarg[ii]);
739  nj = subTempVec.size();
740  for (jj = 0; jj < nj; ++jj)
741  {
742  if (subTempVec[jj] != __notfound)
743  {
744  this->isConstrAct.at(subTempVec[jj]) = true;
745  }
746  }
747  }
748 
749  // Assign Objective Hessian
750 
751  // Assign Constraint Hessian
752 }
RSVScalc.hpp
Provides the infrastructure for calculation of the RSVS equations.
RSVScalctools.hpp
Provides internal tools for RSVS calc.
TriangleActiveDesignVariables
HashedVector< int, int > TriangleActiveDesignVariables(const triangle &triIn, const triangulation &triRSVS, const HashedVector< int, int > &objDvMap, bool useSurfCentreDeriv=true)
Get the active design variables defined by a triangle object.
Definition: RSVScalctools.cpp:43
TrianglePositionDerivatives
void TrianglePositionDerivatives(const triangle &triIn, const triangulation &triRSVS, const HashedVector< int, int > &dvListMap, Eigen::MatrixXd &dPos, Eigen::MatrixXd &HPos, bool useSurfCentreDeriv=true, bool useSurfCentreHessian=false)
Calculate the positional derivatives of a triangle object.
TrianglePointerCoordinates
grid::coordlist TrianglePointerCoordinates(const triangle &triIn, const triangulation &triRSVS)
Get all the coordinates needed to define a triangle.
Definition: RSVScalctools.cpp:18
RSVSmath.hpp
Performs Volume and Area calculations for the RSVS process.
RSVScalc::CalcTriangleEdgeLength
void CalcTriangleEdgeLength(const triangle &triIn, const triangulation &triRSVS, bool isObj=true, bool isConstr=true, bool isDeriv=true)
Calculates the properties of single triangle for 2D RSVS.
Definition: RSVScalc_core.cpp:574
RSVScalc::CalcTriangle
void CalcTriangle(const triangle &triIn, const triangulation &triRSVS, bool isObj=true, bool isConstr=true, bool isDeriv=true)
Calculates the properties of single triangle.
Definition: RSVScalc_core.cpp:32
RSVScalc::CalcTriangleDirectVolume
void CalcTriangleDirectVolume(const triangle &triIn, const triangulation &triRSVS, bool isObj=true, bool isConstr=true, bool isDeriv=true)
Calculates the properties of single triangle using direct calculation.
Definition: RSVScalc_core.cpp:355
RSVScalc::CalcTriangleFD
void CalcTriangleFD(const triangle &triIn, const triangulation &triRSVS, bool isObj=true, bool isConstr=true, bool isDeriv=true)
Calculates the properties of single triangle using Finite difference.
Definition: RSVScalc_core.cpp:199
matrixtools.hpp
Tools to support conversion, display and derivatives of Eigen matrices.
grid::coordlist
std::vector< const std::vector< double > * > coordlist
Defines a list of coordinates.
Definition: mesh.hpp:87
rsvsutils.hpp
Provides various utility functions for the RSVS operation.
snake.hpp
Provides the core restricted surface snake container.
triangulate.hpp
Provides a triangulation for snake and meshes.
vectorarray.hpp
Provides a 2D std::vector based container.
warning.hpp
Provides the error and warning system used by the RSVS3D project.
RSVS3D_ERROR_ARGUMENT
#define RSVS3D_ERROR_ARGUMENT(M)
Throw a invalid_argument.
Definition: warning.hpp:148