CFEASolver CMeshSolver maintenance

Common/include/CConfig.hpp

@@ -992,7 +992,6 @@ class CConfig {
   su2double RefGeom_Penalty,        /*!< \brief Penalty weight value for the reference geometry objective function. */
   RefNode_Penalty,                  /*!< \brief Penalty weight value for the reference node objective function. */
   DV_Penalty;                       /*!< \brief Penalty weight to add a constraint to the total amount of stiffness. */
-  bool addCrossTerm;                /*!< \brief Evaluates the need to add the cross term when setting the adjoint output. */
   unsigned long Nonphys_Points,     /*!< \brief Current number of non-physical points in the solution. */
   Nonphys_Reconstr;                 /*!< \brief Current number of non-physical reconstructions for 2nd-order upwinding. */
   bool ParMETIS;                    /*!< \brief Boolean for activating ParMETIS mode (while testing). */
@@ -2085,7 +2084,7 @@ class CConfig {
    */
   unsigned short GetRefGeom_FileFormat(void) const { return RefGeom_FileFormat; }
 
-    /*!
+  /*!
    * \brief Formulation for 2D elasticity (plane stress - strain)
    * \return Flag to 2D elasticity model.
    */
@@ -2097,27 +2096,16 @@ class CConfig {
    */
   bool GetPrestretch(void) const { return Prestretch; }
 
-  /*!
-    * \brief Decide whether it's necessary to add the cross term for adjoint FSI.
-    * \return <code>TRUE</code> if it's necessary to add the cross term, <code>FALSE</code> otherwise.
-    */
-  bool Add_CrossTerm(void) const { return addCrossTerm; }
-
-  /*!
-    * \brief Set the boolean addCrossTerm to true or false.
-    */
-  void Set_CrossTerm(bool needCrossTerm) { addCrossTerm = needCrossTerm; }
-
   /*!
    * \brief Get the name of the file with the element properties for structural problems.
    * \return Name of the file with the element properties of the structural problem.
    */
   string GetFEA_FileName(void) const { return FEA_FileName; }
 
   /*!
-    * \brief Determine if advanced features are used from the element-based FEA analysis (experimental feature).
-    * \return <code>TRUE</code> is experimental, <code>FALSE</code> is the default behaviour.
-    */
+   * \brief Determine if advanced features are used from the element-based FEA analysis (experimental feature).
+   * \return <code>TRUE</code> is experimental, <code>FALSE</code> is the default behaviour.
+   */
   inline bool GetAdvanced_FEAElementBased(void) const { return FEAAdvancedMode; }
 
   /*!
@@ -8627,7 +8615,7 @@ class CConfig {
    * \param[in] val_index - Index corresponding to the load boundary.
    * \return The pointer to the sine load values.
    */
-  su2double* GetLoad_Sine(void) { return SineLoad_Coeff; }
+  const su2double* GetLoad_Sine(void) const { return SineLoad_Coeff; }
 
   /*!
    * \brief Get the kind of load transfer method we want to use for dynamic problems

Common/include/geometry/elements/CElement.hpp

@@ -449,9 +449,10 @@ class CElement {
    * the latter are needed for compatibility with shape derivatives, there is no problem registering
    * because inactive variables are ignored.
    */
-  inline void SetPreaccIn_Coords(void) {
+  inline void SetPreaccIn_Coords(bool nonlinear = true) {
     AD::SetPreaccIn(RefCoord.data(), nNodes*nDim);
-    AD::SetPreaccIn(CurrentCoord.data(), nNodes*nDim);
+    if (nonlinear)
+      AD::SetPreaccIn(CurrentCoord.data(), nNodes*nDim);
   }
 
   /*!
@@ -462,6 +463,13 @@ class CElement {
     AD::SetPreaccOut(Kt_a.data(), nNodes*nDim);
   }
 
+  /*!
+   * \brief Register the mass matrix as a pre-accumulation output.
+   */
+  inline void SetPreaccOut_Mab(void) {
+    AD::SetPreaccOut(Mab.data(), nNodes*nNodes);
+  }
+
 };
 
 /*!

Common/include/linear_algebra/CSysMatrix.hpp

@@ -705,7 +705,13 @@ class CSysMatrix {
    * \param[in,out] b - The rhs vector (b := b - A_{*,i} * x_i;  b_i = x_i).
    */
   template<class OtherType>
-  void EnforceSolutionAtNode(const unsigned long node_i, const OtherType *x_i, CSysVector<OtherType> & b);
+  void EnforceSolutionAtNode(unsigned long node_i, const OtherType *x_i, CSysVector<OtherType> & b);
+
+  /*!
+   * \brief Version of EnforceSolutionAtNode for a single degree of freedom.
+   */
+  template<class OtherType>
+  void EnforceSolutionAtDOF(unsigned long node_i, unsigned long iVar, OtherType x_i, CSysVector<OtherType> & b);
 
   /*!
    * \brief Sets the diagonal entries of the matrix as the sum of the blocks in the corresponding column.
@@ -718,8 +724,7 @@ class CSysMatrix {
    * \param[in] alpha - The scaling constant.
    * \param[in] B - Matrix being.
    */
-  template<class OtherType>
-  void MatrixMatrixAddition(OtherType alpha, const CSysMatrix<OtherType>& B);
+  void MatrixMatrixAddition(ScalarType alpha, const CSysMatrix& B);
 
   /*!
    * \brief Performs the product of a sparse matrix by a CSysVector.

Common/include/option_structure.hpp

@@ -2020,15 +2020,10 @@ static const MapType<string, ENUM_DIRECTDIFF_VAR> DirectDiff_Var_Map = {
 
 
 enum ENUM_RECORDING {
-  FLOW_CONS_VARS   = 1,
+  SOLUTION_VARIABLES = 1,
   MESH_COORDS = 2,
-  COMBINED    = 3,
-  FEA_DISP_VARS = 4,
-  FLOW_CROSS_TERM = 5,
-  FEM_CROSS_TERM_GEOMETRY = 6,
-  GEOMETRY_CROSS_TERM = 7,
-  ALL_VARIABLES = 8,
-  MESH_DEFORM = 9
+  MESH_DEFORM = 3,
+  SOLUTION_AND_MESH = 4
 };
 
 /*!
@@ -2087,7 +2082,6 @@ enum MPI_QUANTITIES {
   SOLUTION_OLD         =  1,  /*!< \brief Conservative solution old communication. */
   SOLUTION_GRADIENT    =  2,  /*!< \brief Conservative solution gradient communication. */
   SOLUTION_LIMITER     =  3,  /*!< \brief Conservative solution limiter communication. */
-  SOLUTION_DISPONLY    =  4,  /*!< \brief Solution displacement only communication. */
   SOLUTION_PRED        =  5,  /*!< \brief Solution predicted communication. */
   SOLUTION_PRED_OLD    =  6,  /*!< \brief Solution predicted old communication. */
   SOLUTION_GEOMETRY    =  7,  /*!< \brief Geometry solution communication. */

Common/src/interface_interpolation/CNearestNeighbor.cpp

@@ -37,12 +37,7 @@ struct DonorInfo {
   su2double dist;
   unsigned pidx;
   int proc;
-  DonorInfo(su2double d = 0.0, unsigned i = 0, int p = 0) :
-    dist(d), pidx(i), proc(p) { }
-  bool operator< (const DonorInfo& other) const {
-    /*--- Global index is used as tie-breaker to make sorted order independent of initial. ---*/
-    return (dist != other.dist)? (dist < other.dist) : (pidx < other.pidx);
-  }
+  DonorInfo(su2double d = 0.0, unsigned i = 0, int p = 0) : dist(d), pidx(i), proc(p) { }
 };
 
 CNearestNeighbor::CNearestNeighbor(CGeometry ****geometry_container, const CConfig* const* config,  unsigned int iZone,
@@ -139,8 +134,13 @@ void CNearestNeighbor::SetTransferCoeff(const CConfig* const* config) {
         }
       }
 
-      /*--- Find k closest points. ---*/
-      partial_sort(donorInfo.begin(), donorInfo.begin()+nDonor, donorInfo.end());
+      /*--- Find k closest points (need to define the comparator inline or debug build give wrong results). ---*/
+      partial_sort(donorInfo.begin(), donorInfo.begin()+nDonor, donorInfo.end(),
+        [](const DonorInfo& a, const DonorInfo& b) {
+          /*--- Global index is used as tie-breaker to make sorted order independent of initial. ---*/
+          return (a.dist != b.dist)? (a.dist < b.dist) : (a.pidx < b.pidx);
+        }
+      );
 
       /*--- Update stats. ---*/
       numTarget += 1;

Common/src/linear_algebra/CSysMatrix.cpp

@@ -1311,6 +1311,45 @@ void CSysMatrix<ScalarType>::EnforceSolutionAtNode(const unsigned long node_i, c
 
 }
 
+template<class ScalarType>
+template<class OtherType>
+void CSysMatrix<ScalarType>::EnforceSolutionAtDOF(unsigned long node_i, unsigned long iVar,
+                                                  OtherType x_i, CSysVector<OtherType> & b) {
+
+  for (auto index = row_ptr[node_i]; index < row_ptr[node_i+1]; ++index) {
+
+    const auto node_j = col_ind[index];
+
+    /*--- Delete row iVar of block j on row i (bij) and ATTEMPT
+     *    to delete column iVar block i on row j (bji). ---*/
+
+    auto bij = &matrix[index*nVar*nVar];
+    auto bji = GetBlock(node_j, node_i);
+
+    /*--- The "attempt" part. ---*/
+    if (bji != nullptr) {
+      for(auto jVar = 0ul; jVar < nVar; ++jVar) {
+        /*--- Column product. ---*/
+        b[node_j*nVar+jVar] -= bji[jVar*nVar+iVar] * x_i;
+        /*--- Delete entries. ---*/
+        bji[jVar*nVar+iVar] = 0.0;
+      }
+    }
+
+    /*--- Delete row. ---*/
+    for(auto jVar = 0ul; jVar < nVar; ++jVar)
+      bij[iVar*nVar+jVar] = 0.0;
+
+    /*--- Set the diagonal entry of the block to 1. ---*/
+    if (node_j == node_i)
+      bij[iVar*(nVar+1)] = 1.0;
+  }
+
+  /*--- Set know solution in rhs vector. ---*/
+  b(node_i, iVar) = x_i;
+
+}
+
 template<class ScalarType>
 void CSysMatrix<ScalarType>::SetDiagonalAsColumnSum() {
 
@@ -1329,8 +1368,7 @@ void CSysMatrix<ScalarType>::SetDiagonalAsColumnSum() {
 }
 
 template<class ScalarType>
-template<class OtherType>
-void CSysMatrix<ScalarType>::MatrixMatrixAddition(OtherType alpha, const CSysMatrix<OtherType>& B) {
+void CSysMatrix<ScalarType>::MatrixMatrixAddition(ScalarType alpha, const CSysMatrix<ScalarType>& B) {
 
   /*--- Check that the sparse structure is shared between the two matrices,
    *    comparing pointers is ok as they are obtained from CGeometry. ---*/
@@ -1344,7 +1382,7 @@ void CSysMatrix<ScalarType>::MatrixMatrixAddition(OtherType alpha, const CSysMat
 
   SU2_OMP_FOR_STAT(omp_light_size)
   for (auto i = 0ul; i < nnz*nVar*nEqn; ++i)
-    matrix[i] += PassiveAssign<ScalarType,OtherType>(alpha*B.matrix[i]);
+    matrix[i] += alpha*B.matrix[i];
 
 }
 
@@ -1403,7 +1441,7 @@ template class CSysMatrix<su2double>;
 template void  CSysMatrix<su2double>::InitiateComms(const CSysVector<su2double>&, CGeometry*, CConfig*, unsigned short) const;
 template void  CSysMatrix<su2double>::CompleteComms(CSysVector<su2double>&, CGeometry*, CConfig*, unsigned short) const;
 template void  CSysMatrix<su2double>::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector<su2double>&);
-template void  CSysMatrix<su2double>::MatrixMatrixAddition(su2double, const CSysMatrix<su2double>&);
+template void  CSysMatrix<su2double>::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double, CSysVector<su2double>&);
 
 #ifdef CODI_REVERSE_TYPE
 template class CSysMatrix<passivedouble>;
@@ -1413,6 +1451,6 @@ template void  CSysMatrix<passivedouble>::CompleteComms(CSysVector<passivedouble
 template void  CSysMatrix<passivedouble>::CompleteComms(CSysVector<su2double>&, CGeometry*, CConfig*, unsigned short) const;
 template void  CSysMatrix<passivedouble>::EnforceSolutionAtNode(unsigned long, const passivedouble*, CSysVector<passivedouble>&);
 template void  CSysMatrix<passivedouble>::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector<su2double>&);
-template void  CSysMatrix<passivedouble>::MatrixMatrixAddition(passivedouble, const CSysMatrix<passivedouble>&);
-template void  CSysMatrix<passivedouble>::MatrixMatrixAddition(su2double, const CSysMatrix<su2double>&);
+template void  CSysMatrix<passivedouble>::EnforceSolutionAtDOF(unsigned long, unsigned long, passivedouble, CSysVector<passivedouble>&);
+template void  CSysMatrix<passivedouble>::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double, CSysVector<su2double>&);
 #endif

SU2_CFD/include/iteration_structure.hpp

@@ -296,12 +296,6 @@ class CIteration {
                                  unsigned short iZone,
                                  unsigned short iInst){}
 
-  virtual void InitializeAdjoint_CrossTerm(CSolver *****solver,
-                                 CGeometry ****geometry,
-                                 CConfig **config,
-                                 unsigned short iZone,
-                                 unsigned short iInst){}
-
   virtual void RegisterInput(CSolver *****solver,
                              CGeometry ****geometry,
                              CConfig** config,
@@ -1195,7 +1189,7 @@ class CDiscAdjFluidIteration : public CIteration {
    * \param[in] config - Definition of the particular problem.
    * \param[in] iZone - Index of the zone.
    * \param[in] iInst - Index of the instance.
-   * \param[in] kind_recording - Kind of recording, either FLOW_CONS_VARS or MESH_COORDS
+   * \param[in] kind_recording - Kind of recording, either FLOW_SOLUTION_VARIABLES or MESH_COORDS
    */
   void RegisterInput(CSolver *****solver,
                      CGeometry ****geometry,
@@ -1219,16 +1213,6 @@ class CDiscAdjFluidIteration : public CIteration {
                       unsigned short iZone,
                       unsigned short iInst);
 
-  /*!
-   * \brief Initializes the adjoints of the output variables of the meanflow iteration - without the contribution of the objective function
-   * \param[in] solver - Container vector with all the solutions.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] iZone - Index of the zone.
-   * \param[in] iInst - Index of the instance.
-   */
-  void InitializeAdjoint_CrossTerm(CSolver *****solver, CGeometry ****geometry, CConfig **config, unsigned short iZone, unsigned short iInst);
-
   /*!
    * \brief Record a single iteration of the direct mean flow system.
    * \param[in] output - Pointer to the COutput class.
@@ -1499,16 +1483,6 @@ class CDiscAdjFEAIteration : public CIteration {
    */
   void InitializeAdjoint(CSolver *****solver, CGeometry ****geometry, CConfig** config, unsigned short iZone, unsigned short iInst);
 
-  /*!
-   * \brief Initializes the adjoints of the output variables of the FEM iteration - without the contribution of the objective function
-   * \param[in] solver - Container vector with all the solutions.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] iZone - Index of the zone.
-   * \param[in] iInst - Index of the zone.
-   */
-  void InitializeAdjoint_CrossTerm(CSolver *****solver, CGeometry ****geometry, CConfig **config, unsigned short iZone, unsigned short iInst);
-
   /*!
    * \brief Record a single iteration of the direct FEM system.
    * \param[in] output - Pointer to the COutput class.