Discrete adjoint for deforming meshes

Common/include/geometry/dual_grid/CVertex.hpp

@@ -37,26 +37,23 @@
  */
 class CVertex : public CDualGrid {
 protected:
-  unsigned long *Nodes; /*!< \brief Vector to store the global nodes of an element. */
-  su2double *Normal;      /*!< \brief Normal coordinates of the element and its center of gravity. */
-  su2double Aux_Var;      /*!< \brief Auxiliar variable defined only on the surface. */
-  su2double CartCoord[3];   /*!< \brief Vertex cartesians coordinates. */
-  su2double VarCoord[3];    /*!< \brief Used for storing the coordinate variation due to a surface modification. */
-  su2double *VarRot;   /*!< \brief Used for storing the rotation variation due to a surface modification. */
-  long PeriodicPoint[5];      /*!< \brief Store the periodic point of a boundary (iProcessor, iPoint) */
-  bool ActDisk_Perimeter;     /*!< \brief Identify nodes at the perimeter of the actuator disk */
-  short Rotation_Type;      /*!< \brief Type of rotation associated with the vertex (MPI and periodic) */
+  unsigned long Nodes[1];        /*!< \brief Vector to store the global nodes of an element. */
+  su2double Normal[3];           /*!< \brief Normal coordinates of the element and its center of gravity. */
+  su2double Aux_Var;             /*!< \brief Auxiliar variable defined only on the surface. */
+  su2double CartCoord[3];        /*!< \brief Vertex cartesians coordinates. */
+  su2double VarCoord[3];         /*!< \brief Used for storing the coordinate variation due to a surface modification. */
+  su2double *VarRot;             /*!< \brief Used for storing the rotation variation due to a surface modification. */
+  long PeriodicPoint[5];         /*!< \brief Store the periodic point of a boundary (iProcessor, iPoint) */
+  bool ActDisk_Perimeter;        /*!< \brief Identify nodes at the perimeter of the actuator disk */
+  short Rotation_Type;           /*!< \brief Type of rotation associated with the vertex (MPI and periodic) */
   unsigned long Normal_Neighbor; /*!< \brief Index of the closest neighbor. */
-  unsigned long *Donor_Points; /*!< \brief indices of donor points for interpolation across zones */
-  unsigned long *Donor_Proc; /*!< \brief indices of donor processor for interpolation across zones in parallel */
-  unsigned long Donor_Elem;   /*!< \brief Store the donor element for interpolation across zones/ */
-  unsigned short Donor_Face;  /*!<\brief Store the donor face (w/in donor element) for interpolation across zones */
-  su2double Basis_Function[3]; /*!< \brief Basis function values for interpolation across zones. */
-  su2double *Donor_Coeff; /*!\brief Store a list of coefficients corresponding to the donor points. */
-  unsigned short nDonor_Points; /*!\brief Number of points in Donor_Points; at least there will be one donor point (if the mesh is matching)*/
+  unsigned long *Donor_Points;   /*!< \brief indices of donor points for interpolation across zones */
+  unsigned long *Donor_Proc;     /*!< \brief indices of donor processor for interpolation across zones in parallel */
+  su2double Basis_Function[3];   /*!< \brief Basis function values for interpolation across zones. */
+  su2double *Donor_Coeff;        /*!< \brief Store a list of coefficients corresponding to the donor points. */
+  unsigned short nDonor_Points;  /*!< \brief Number of points in Donor_Coeff. */
 
 public:
-
   /*!
    * \brief Constructor of the class.
    * \param[in] val_point - Node of the vertex.
@@ -113,6 +110,11 @@ class CVertex : public CDualGrid {
    */
   inline su2double *GetNormal(void) override { return Normal; }
 
+  /*!
+   * \brief Get the ith component of the normal.
+   */
+  inline su2double GetNormal(unsigned short iDim) const { return Normal[iDim]; }
+
   /*!
    * \brief Initialize normal vector.
    */
@@ -308,30 +310,6 @@ class CVertex : public CDualGrid {
    */
   inline bool GetActDisk_Perimeter(void) const { return ActDisk_Perimeter; }
 
-  /*!
-   * \brief Set the donor element of a vertex for interpolation across zones.
-   * \param[in] val_donorelem - donor element index.
-   */
-  inline void SetDonorElem(long val_donorelem) { Donor_Elem = val_donorelem; }
-
-  /*!
-   * \brief Get the donor element of a vertex for interpolation across zones.
-   * \return Value of the donor element of a vertex.
-   */
-  inline long GetDonorElem(void) const { return Donor_Elem; }
-
-  /*!
-   * \brief Set the donor face of a vertex for interpolation across zones.
-   * \param[in] val_donorface- donor face index (w/in donor elem).
-   */
-  inline void SetDonorFace(unsigned short val_donorface) { Donor_Face = val_donorface; }
-
-  /*!
-   * \brief Get the donor face of a vertex for interpolation across zones.
-   * \return Value of the donor face index (w/in donor elem).
-   */
-  inline unsigned short GetDonorFace(void) const { return Donor_Face; }
-
   /*!
    * \brief Set the finite element basis functions needed for interpolation.
    * \param[in] val_node - a node index of the owner element.
@@ -358,17 +336,6 @@ class CVertex : public CDualGrid {
    */
   inline unsigned long GetNormal_Neighbor(void) const { return Normal_Neighbor; }
 
-  /*!
-   * \brief Increment the number of donor points by 1.
-   */
-  inline void IncrementnDonor(void) {nDonor_Points++;}
-
-  /*!
-   * \brief Set the value of nDonor_Points
-   * \param[in] nDonor - the number of donor points
-   */
-  inline void SetnDonorPoints(unsigned short nDonor) {nDonor_Points = nDonor;}
-
   /*!
    * \brief Return the value of nDonor_Points
    * \return nDonor - the number of donor points
@@ -419,9 +386,9 @@ class CVertex : public CDualGrid {
 
   /*!
    * \brief Allocate memory based on how many donor points need to be stored.
-   * Uses nDonor_Points
+   * \param[in] nDonor - the number of donor points
    */
-  void Allocate_DonorInfo(void);
+  void Allocate_DonorInfo(unsigned short nDonor);
 
   /*!
    * \brief Get the rotation variation

Common/include/geometry/elements/CElement.hpp

@@ -7,7 +7,7 @@
  *
  * SU2 Project Website: https://su2code.github.io
  *
- * The SU2 Project is maintained by the SU2 Foundation 
+ * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
  * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
@@ -627,6 +627,26 @@ class CQUAD4 final : public CElementWithKnownSizes<4,4,2> {
    */
   CQUAD4();
 
+  /*!
+   * \brief Shape functions (Ni) evaluated at point Xi,Eta.
+   */
+  inline static void ShapeFunctions(su2double Xi, su2double Eta, su2double* Ni) {
+    Ni[0] = 0.25*(1.0-Xi)*(1.0-Eta);
+    Ni[1] = 0.25*(1.0+Xi)*(1.0-Eta);
+    Ni[2] = 0.25*(1.0+Xi)*(1.0+Eta);
+    Ni[3] = 0.25*(1.0-Xi)*(1.0+Eta);
+  }
+
+  /*!
+   * \brief Shape function Jacobian (dNi) evaluated at point Xi,Eta.
+   */
+  inline static void ShapeFunctionJacobian(su2double Xi, su2double Eta, su2double dNi[][2]) {
+    dNi[0][0] = -0.25*(1.0-Eta);  dNi[0][1] = -0.25*(1.0-Xi);
+    dNi[1][0] =  0.25*(1.0-Eta);  dNi[1][1] = -0.25*(1.0+Xi);
+    dNi[2][0] =  0.25*(1.0+Eta);  dNi[2][1] =  0.25*(1.0+Xi);
+    dNi[3][0] = -0.25*(1.0+Eta);  dNi[3][1] =  0.25*(1.0-Xi);
+  }
+
   /*!
    * \brief Compute the value of the area of the element.
    * \param[in] mode - Type of coordinates to consider in the computation.

Common/include/grid_movement_structure.hpp

@@ -1290,183 +1290,6 @@ class CVolumetricMovement : public CGridMovement {
   virtual void Boundary_Dependencies(CGeometry **geometry, CConfig *config);
 };
 
-/*! 
- * \class CElasticityMovement
- * \brief Class for moving the volumetric numerical grid using the new linear elasticity solver.
- * \author R.Sanchez, based on CVolumetricMovement developments of F. Palacios, A. Bueno, T. Economon, S. Padron
- * \version 7.0.3 "Blackbird"
- */
-class CElasticityMovement : public CVolumetricMovement {
-protected:
-
-  unsigned short nDim;    /*!< \brief Number of dimensions. */
-  unsigned short nVar;    /*!< \brief Number of variables. */
-
-  unsigned long nPoint;   /*!< \brief Number of points. */
-  unsigned long nPointDomain;   /*!< \brief Number of points in the domain. */
-
-  unsigned long nIterMesh;   /*!< \brief Number of iterations in the mesh update. +*/
-  su2double valResidual;
-
-  su2double *Residual,         /*!< \brief Auxiliary nDim vector. */
-  *Solution;         /*!< \brief Auxiliary nDim vector. */
-
-  su2double **matrixZeros;     /*!< \brief Submatrix to make zeros and impose boundary conditions. */
-  su2double **matrixId;        /*!< \brief Diagonal submatrix to impose boundary conditions. */
-
-  su2double MinVolume;
-  su2double MaxVolume;
-
-#ifndef CODI_FORWARD_TYPE
-  CSysSolve<passivedouble>  System;
-  CSysMatrix<passivedouble> StiffMatrix; /*!< \brief Matrix to store the point-to-point stiffness. */
-#else
-  CSysSolve<su2double>  System;
-  CSysMatrix<su2double> StiffMatrix;
-#endif
-  CSysVector<su2double> LinSysSol;
-  CSysVector<su2double> LinSysRes;
-
-  su2double E;                  /*!< \brief Young's modulus of elasticity. */
-  su2double Nu;                 /*!< \brief Poisson's ratio. */
-
-  su2double Mu;                 /*!< \brief Lame's coeficient. */
-  su2double Lambda;             /*!< \brief Lame's coeficient. */
-
-  su2double **Jacobian_ij;      /*!< \brief Submatrix to store the constitutive term for node ij. */
-
-  su2double **Ba_Mat,           /*!< \brief Matrix B for node a - Auxiliary. */
-  **Bb_Mat;           /*!< \brief Matrix B for node b - Auxiliary. */
-  su2double **KAux_ab;          /*!< \brief Stiffness sub-term  - Auxiliary. */
-  su2double **D_Mat;            /*!< \brief Constitutive matrix - Auxiliary. */
-  su2double **GradNi_Ref_Mat;   /*!< \brief Gradients of Ni - Auxiliary. */
-
-public:
-
-  CElement** element_container;  /*!< \brief Container which stores the element information. */
-
-  /*!
-   * \brief Constructor of the class.
-   */
-  CElasticityMovement(CGeometry *geometry, CConfig *config);
-
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CElasticityMovement(void);
-
-  /*!
-   * \brief Grid deformation using the linear elasticity equations.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] UpdateGeo - Update geometry.
-   * \param[in] Derivative - Compute the derivative (disabled by default). Does not actually deform the grid if enabled.
-   */
-  void SetVolume_Deformation_Elas(CGeometry *geometry, CConfig *config, bool UpdateGeo, bool screen_output, bool Derivative = false);
-
-  /*!
-   * \brief Update the value of the coordinates after the grid movement.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void UpdateGridCoord(CGeometry *geometry, CConfig *config);
-
-  /*!
-   * \brief Update the dual grid after the grid movement (edges and control volumes).
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void UpdateDualGrid(CGeometry *geometry, CConfig *config);
-
-  /*!
-   * \brief Update the coarse multigrid levels after the grid movement.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void UpdateMultiGrid(CGeometry **geometry, CConfig *config);
-
-  /*!
-   * \brief Check the boundary vertex that are going to be moved.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void SetBoundaryDisplacements(CGeometry *geometry, CConfig *config);
-
-  /*!
-   * \brief Set the boundary displacements to 0.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] val_marker -
-   */
-  void SetClamped_Boundary(CGeometry *geometry, CConfig *config, unsigned short val_marker);
-
-  /*!
-   * \brief Set the boundary displacements to the imposed external value.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void SetMoving_Boundary(CGeometry *geometry, CConfig *config, unsigned short val_marker);
-
-  /*!
-   * \brief Compute the min and max volume for the stiffness matrix for grid deformation.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \return Value of the length of the smallest edge of the grid.
-   */
-  void SetMinMaxVolume(CGeometry *geometry, CConfig *config);
-
-  /*!
-   * \brief Compute the min and max volume for the stiffness matrix for grid deformation.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \return Value of the length of the smallest edge of the grid.
-   */
-  void SetStiffnessMatrix(CGeometry *geometry, CConfig *config);
-
-  /*!
-   * \brief Store the number of iterations when moving the mesh.
-   * \param[in] val_nIterMesh - Number of iterations.
-   */
-  void Set_nIterMesh(unsigned long val_nIterMesh);
-
-  /*!
-   * \brief Retrieve the number of iterations when moving the mesh.
-   * \param[out] Number of iterations.
-   */
-  unsigned long Get_nIterMesh(void);
-
-  /*!
-   * \brief Compute the stiffness of the element and the parameters Lambda and Mu
-   */
-  void Set_Element_Stiffness(su2double ElemVolume, CConfig *config);
-
-  /*!
-   * \brief Compute the stiffness of the element and the parameters Lambda and Mu
-   */
-  void Compute_Element_Contribution(CElement *element, CConfig *config);
-
-  /*!
-   * \brief Compute the constitutive matrix in an element for mesh deformation problems
-   * \param[in] element_container - Element structure for the particular element integrated.
-   */
-  void Compute_Constitutive_Matrix(void);
-
-  /*!
-   * \brief Set the boundary displacements in the mesh side of the problem
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void Transfer_Boundary_Displacements(CGeometry *geometry, CConfig *config, unsigned short val_marker);
-
-  /*!
-   * \brief Set the boundary displacements in the mesh side of the problem
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void Boundary_Dependencies(CGeometry **geometry, CConfig *config);
-
-};
-
 /*!
  * \class CSurfaceMovement
  * \brief Class for moving the surface numerical grid.
@@ -1571,46 +1394,6 @@ class CSurfaceMovement : public CGridMovement {
    */
   void Moving_Walls(CGeometry *geometry, CConfig *config, unsigned short iZone, unsigned long iter);
 
-  /*!
-   * \brief Computes the displacement of a translating surface for a dynamic mesh simulation.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] iter - Current physical time iteration.
-   * \param[in] iZone - Zone number in the mesh.
-   */
-  void Surface_Translating(CGeometry *geometry, CConfig *config,
-                           unsigned long iter, unsigned short iZone);
-
-  /*!
-   * \brief Computes the displacement of a plunging surface for a dynamic mesh simulation.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] iter - Current physical time iteration.
-   * \param[in] iZone - Zone number in the mesh.
-   */
-  void Surface_Plunging(CGeometry *geometry, CConfig *config,
-                        unsigned long iter, unsigned short iZone);
-
-  /*!
-   * \brief Computes the displacement of a pitching surface for a dynamic mesh simulation.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] iter - Current physical time iteration.
-   * \param[in] iZone - Zone number in the mesh.
-   */
-  void Surface_Pitching(CGeometry *geometry, CConfig *config,
-                        unsigned long iter, unsigned short iZone);
-
-  /*!
-   * \brief Computes the displacement of a rotating surface for a dynamic mesh simulation.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] iter - Current physical time iteration.
-   * \param[in] iZone - Zone number in the mesh.
-   */
-  void Surface_Rotating(CGeometry *geometry, CConfig *config,
-                        unsigned long iter, unsigned short iZone);
-
   /*!
    * \brief Computes the displacement of a rotating surface for a dynamic mesh simulation.
    * \param[in] geometry - Geometrical definition of the problem.

Common/include/grid_movement_structure.inl

@@ -216,10 +216,6 @@ inline void CVolumetricMovement::SetVolume_Deformation_Elas(CGeometry *geometry,
 
 inline void CVolumetricMovement::Boundary_Dependencies(CGeometry **geometry, CConfig *config) {  }
 
-inline void CElasticityMovement::Set_nIterMesh(unsigned long val_nIterMesh) { nIterMesh = val_nIterMesh; }
-
-inline unsigned long CElasticityMovement::Get_nIterMesh() { return nIterMesh; }
-
 inline bool CSurfaceMovement::CheckFFDBoxDefinition(CConfig *config, unsigned short iDV) {
   for (unsigned short iFFDBox = 0; iFFDBox < GetnFFDBox(); iFFDBox++) {
     if (FFDBox[iFFDBox]->GetTag() == config->GetFFDTag(iDV)) { return true;}