diff --git a/.gitignore b/.gitignore index 2177d44e12e1..de46aa105d4c 100644 --- a/.gitignore +++ b/.gitignore @@ -103,3 +103,10 @@ su2preconfig.timestamp # Clangd server files .cache + + + +ninja-win.zip +ninja.exe +.gitignore +Docs/html diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index c28c3aaef54e..8670796a3db2 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -189,6 +189,7 @@ class CConfig { nMarker_ActDiskBemOutlet_Axis, /*!< \brief Number of actuator disk BEM outlet markers passed to MARKER_ACTDISK_BEM_AXIS. */ nMarker_Deform_Mesh_Sym_Plane, /*!< \brief Number of markers with symmetric deformation */ nMarker_Deform_Mesh, /*!< \brief Number of deformable markers at the boundary. */ + nMarker_Deform_Mesh_Internal, /*!< \brief Number of internal markers allowed to freely deform. */ nMarker_Fluid_Load, /*!< \brief Number of markers in which the flow load is computed/employed. */ nMarker_Fluid_InterfaceBound, /*!< \brief Number of fluid interface markers. */ nMarker_CHTInterface, /*!< \brief Number of conjugate heat transfer interface markers. */ @@ -238,6 +239,7 @@ class CConfig { *Marker_NearFieldBound, /*!< \brief Near Field boundaries markers. */ *Marker_Deform_Mesh, /*!< \brief Deformable markers at the boundary. */ *Marker_Deform_Mesh_Sym_Plane, /*!< \brief Marker with symmetric deformation. */ + *Marker_Deform_Mesh_Internal, /*!< \brief Internal marker allowed to freely deform. */ *Marker_Fluid_Load, /*!< \brief Markers in which the flow load is computed/employed. */ *Marker_Fluid_InterfaceBound, /*!< \brief Fluid interface markers. */ *Marker_CHTInterface, /*!< \brief Conjugate heat transfer interface markers. */ @@ -652,6 +654,10 @@ class CConfig { su2double Deform_Tol_Factor; /*!< \brief Factor to multiply smallest volume for deform tolerance (0.001 default) */ su2double Deform_Coeff; /*!< \brief Deform coeffienct */ su2double Deform_Limit; /*!< \brief Deform limit */ + DEFORM_KIND Deform_Kind; /*!< \brief Type of mesh deformation */ + bool RBF_DataReduction; /*!< \brief Determines use of data reduction methods for RBF mesh deformation. */ + su2double RBF_GreedyTolerance; /*!< \brief Tolerance used in the greedy data reduction for RBF mesh deformation. */ + su2double RBF_GreedyCorrectionFactor; /*!< \brief Correction factor used in the greedy algorithm for RBF mesh deformation. */ unsigned short FFD_Continuity; /*!< \brief Surface continuity at the intersection with the FFD */ unsigned short FFD_CoordSystem; /*!< \brief Define the coordinates system */ su2double Deform_ElasticityMod, /*!< \brief Young's modulus for volume deformation stiffness model */ @@ -746,6 +752,7 @@ class CConfig { *Marker_All_Moving, /*!< \brief Global index for moving surfaces using the grid information. */ *Marker_All_Deform_Mesh, /*!< \brief Global index for deformable markers at the boundary. */ *Marker_All_Deform_Mesh_Sym_Plane, /*!< \brief Global index for markers with symmetric deformations. */ + *Marker_All_Deform_Mesh_Internal, /*!< \brief Global index for internal markers with free deformation. */ *Marker_All_Fluid_Load, /*!< \brief Global index for markers in which the flow load is computed/employed. */ *Marker_All_PyCustom, /*!< \brief Global index for Python customizable surfaces using the grid information. */ *Marker_All_Designing, /*!< \brief Global index for moving using the grid information. */ @@ -762,6 +769,7 @@ class CConfig { *Marker_CfgFile_Moving, /*!< \brief Global index for moving surfaces using the config information. */ *Marker_CfgFile_Deform_Mesh, /*!< \brief Global index for deformable markers at the boundary. */ *Marker_CfgFile_Deform_Mesh_Sym_Plane, /*!< \brief Global index for markers with symmetric deformations. */ + *Marker_CfgFile_Deform_Mesh_Internal, /*!< \brief Global index for internal markers with free deformation. */ *Marker_CfgFile_Fluid_Load, /*!< \brief Global index for markers in which the flow load is computed/employed. */ *Marker_CfgFile_PyCustom, /*!< \brief Global index for Python customizable surfaces using the config information. */ *Marker_CfgFile_DV, /*!< \brief Global index for design variable markers using the config information. */ @@ -3487,6 +3495,14 @@ class CConfig { */ void SetMarker_All_Deform_Mesh_Sym_Plane(unsigned short val_marker, unsigned short val_deform) { Marker_All_Deform_Mesh_Sym_Plane[val_marker] = val_deform; } + /*! + * \brief Set if a marker val_marker allows deformation at the boundary. + * \param[in] val_marker - Index of the marker in which we are interested. + * \param[in] val_interface - 0 or 1 depending if the the marker is or not a DEFORM_MESH_SYM_PLANE marker. + */ + void SetMarker_All_Deform_Mesh_Internal(unsigned short val_marker, unsigned short val_deform) { Marker_All_Deform_Mesh_Internal[val_marker] = val_deform; } + + /*! * \brief Set if a in marker val_marker the flow load will be computed/employed. * \param[in] val_marker - Index of the marker in which we are interested. @@ -3637,6 +3653,13 @@ class CConfig { * \return 0 or 1 depending if the marker belongs to the DEFORM_MESH_SYM_PLANE subset. */ unsigned short GetMarker_All_Deform_Mesh_Sym_Plane(unsigned short val_marker) const { return Marker_All_Deform_Mesh_Sym_Plane[val_marker]; } + + /*! + * \brief Get whether marker val_marker is a DEFORM_MESH_SYM_PLANE marker + * \param[in] val_marker - 0 or 1 depending if the the marker belongs to the DEFORM_MESH_SYM_PLANE subset. + * \return 0 or 1 depending if the marker belongs to the DEFORM_MESH_SYM_PLANE subset. + */ + unsigned short GetMarker_All_Deform_Mesh_Internal(unsigned short val_marker) const { return Marker_All_Deform_Mesh_Internal[val_marker]; } /*! * \brief Get whether marker val_marker is a Fluid_Load marker @@ -4375,6 +4398,30 @@ class CConfig { */ bool GetFFD_Symmetry_Plane(void) const { return FFD_Symmetry_Plane; } + /*! + * \brief Get the type of mesh deformation method. + * \return type of mesh deformation. + */ + DEFORM_KIND GetDeform_Kind() const { return Deform_Kind; } + + /*! + * \brief Determines use of data reduction methods for RBF mesh deformation. + * \return TRUE means that data reduction is used. + */ + bool GetRBF_DataReduction(void) const { return RBF_DataReduction; } + + /*! + * \brief Determines use of data reduction methods for RBF mesh deformation. + * \return TRUE means that data reduction is used. + */ + su2double GetRBF_DataRedTolerance(void) const { return RBF_GreedyTolerance; } + + /*! + * \brief Determines use of data reduction methods for RBF mesh deformation. + * \return TRUE means that data reduction is used. + */ + su2double GetRBF_DataRedCorrectionFactor(void) const { return RBF_GreedyCorrectionFactor; } + /*! * \brief Get the kind of SU2 software component. * \return Kind of the SU2 software component. @@ -6359,6 +6406,12 @@ class CConfig { */ unsigned short GetMarker_CfgFile_Deform_Mesh_Sym_Plane(const string& val_marker) const; + /*! + * \brief Get the DEFORM_MESH_INTERNAL information from the config definition for the marker val_marker. + * \return DEFORM_MESH_INTERNAL information of the boundary in the config information for the marker val_marker. + */ + unsigned short GetMarker_CfgFile_Deform_Mesh_Internal(const string& val_marker) const; + /*! * \brief Get the Fluid_Load information from the config definition for the marker val_marker. * \return Fluid_Load information of the boundary in the config information for the marker val_marker. @@ -6714,6 +6767,12 @@ class CConfig { */ unsigned short GetMarker_Deform_Mesh_Sym_Plane(const string& val_marker) const; + /*! + * \brief Get the internal index for a DEFORM_MESH_SYM_PLANE boundary val_marker. + * \return Internal index for a DEFORM_MESH_SYM_PLANE boundary val_marker. + */ + unsigned short GetMarker_Deform_Mesh_Internal(const string& val_marker) const; + /*! * \brief Get a bool for whether the marker is deformed. val_marker. * \param[in] val_marker - Name of the marker to test. diff --git a/Common/include/grid_movement/CLinearElasticity.hpp b/Common/include/grid_movement/CLinearElasticity.hpp new file mode 100644 index 000000000000..e1745b06db28 --- /dev/null +++ b/Common/include/grid_movement/CLinearElasticity.hpp @@ -0,0 +1,246 @@ +/*! + * \file CLinearElasticity.hpp + * \brief Headers of the CLinearElasticity class. + * \author F. Palacios, A. Bueno, T. Economon, S. Padron. + * \version 8.0.1 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once +#include "CVolumetricMovement.hpp" +#include "../linear_algebra/CSysMatrix.hpp" +#include "../linear_algebra/CSysVector.hpp" +#include "../linear_algebra/CSysSolve.hpp" + +/*! + * \class CLinearElasticity + * \brief Class for moving the volumetric numerical grid using the linear elasticity analogy. + * \author F. Palacios, A. Bueno, T. Economon, S. Padron. + */ + +class CLinearElasticity final: public CVolumetricMovement{ + protected: + 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. +*/ + + + #ifndef CODI_FORWARD_TYPE + CSysMatrix StiffMatrix; /*!< \brief Stiffness matrix of the elasticity problem. */ + CSysSolve System; /*!< \brief Linear solver/smoother. */ +#else + CSysMatrix StiffMatrix; + CSysSolve System; +#endif + CSysVector LinSysSol; + CSysVector LinSysRes; + + public: + /*! + * \brief Constructor of the class. + */ + CLinearElasticity(CGeometry* geometry, CConfig* config); + + /*! + * \brief Destructor of the class. + */ + ~CLinearElasticity() override; + + /*! + * \brief Grid deformation using the spring analogy method. + * \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(CGeometry* geometry, CConfig* config, bool UpdateGeo, bool Derivative, + bool ForwardProjectionDerivative); + + /*! + * \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 derivatives 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_Derivatives(CGeometry* geometry, CConfig* config, bool ForwardProjectionDerivative); + + /*! + * \brief Compute the minimum distance to the nearest solid surface. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void ComputeSolid_Wall_Distance(CGeometry* geometry, CConfig* config, su2double& MinDistance, + su2double& MaxDistance) const; + + /*! + * \brief Compute the stiffness matrix for grid deformation using spring analogy. + * \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. + */ + su2double SetFEAMethodContributions_Elem(CGeometry* geometry, CConfig* config); + + /*! + * \brief Build the stiffness matrix for a 3-D hexahedron element. The result will be placed in StiffMatrix_Elem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] StiffMatrix_Elem - Element stiffness matrix to be filled. + * \param[in] CoordCorners - Index value for Node 1 of the current hexahedron. + * \param[in] PointCorners - Index values for element corners + * \param[in] nNodes - Number of nodes defining the element. + * \param[in] scale + */ + void SetFEA_StiffMatrix2D(CGeometry* geometry, CConfig* config, su2double** StiffMatrix_Elem, + unsigned long PointCorners[8], su2double CoordCorners[8][3], + unsigned short nNodes, su2double ElemVolume, su2double ElemDistance) ; + + /*! + * \brief Build the stiffness matrix for a 3-D hexahedron element. The result will be placed in StiffMatrix_Elem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] StiffMatrix_Elem - Element stiffness matrix to be filled. + * \param[in] CoordCorners - Index value for Node 1 of the current hexahedron. + * \param[in] PointCorners - Index values for element corners + * \param[in] nNodes - Number of nodes defining the element. + * \param[in] scale + */ + void SetFEA_StiffMatrix3D(CGeometry* geometry, CConfig* config, su2double** StiffMatrix_Elem, + unsigned long PointCorners[8], su2double CoordCorners[8][3], + unsigned short nNodes, su2double ElemVolume, su2double ElemDistance); + + /*! + * \brief Add the stiffness matrix for a 2-D triangular element to the global stiffness matrix for the entire mesh + * (node-based). \param[in] geometry - Geometrical definition of the problem. \param[in] StiffMatrix_Elem - Element + * stiffness matrix to be filled. \param[in] PointCorners - Index values for element corners \param[in] nNodes - + * Number of nodes defining the element. + */ + void AddFEA_StiffMatrix(CGeometry* geometry, su2double** StiffMatrix_Elem, + unsigned long PointCorners[8], unsigned short nNodes); + + /*! + * \brief Shape functions and derivative of the shape functions + * \param[in] Xi - Local coordinates. + * \param[in] Eta - Local coordinates. + * \param[in] Zeta - Local coordinates. + * \param[in] CoordCorners - Coordiantes of the corners. + * \param[in] DShapeFunction - Shape function information + */ + su2double ShapeFunc_Hexa(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3], + su2double DShapeFunction[8][4]); + + /*! + * \brief Shape functions and derivative of the shape functions + * \param[in] Xi - Local coordinates. + * \param[in] Eta - Local coordinates. + * \param[in] Zeta - Local coordinates. + * \param[in] CoordCorners - Coordiantes of the corners. + * \param[in] DShapeFunction - Shape function information + */ + su2double ShapeFunc_Tetra(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3], + su2double DShapeFunction[8][4]); + + /*! + * \brief Shape functions and derivative of the shape functions + * \param[in] Xi - Local coordinates. + * \param[in] Eta - Local coordinates. + * \param[in] Zeta - Local coordinates. + * \param[in] CoordCorners - Coordiantes of the corners. + * \param[in] DShapeFunction - Shape function information + */ + su2double ShapeFunc_Pyram(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3], + su2double DShapeFunction[8][4]); + + /*! + * \brief Shape functions and derivative of the shape functions + * \param[in] Xi - Local coordinates. + * \param[in] Eta - Local coordinates. + * \param[in] Zeta - Local coordinates. + * \param[in] CoordCorners - Coordiantes of the corners. + * \param[in] DShapeFunction - Shape function information + */ + su2double ShapeFunc_Prism(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3], + su2double DShapeFunction[8][4]); + + /*! + * \brief Shape functions and derivative of the shape functions + * \param[in] Xi - Local coordinates. + * \param[in] Eta - Local coordinates. + * \param[in] CoordCorners - Coordiantes of the corners. + * \param[in] DShapeFunction - Shape function information + */ + su2double ShapeFunc_Triangle(su2double Xi, su2double Eta, su2double CoordCorners[8][3], + su2double DShapeFunction[8][4]); + + /*! + * \brief Shape functions and derivative of the shape functions + * \param[in] Xi - Local coordinates. + * \param[in] Eta - Local coordinates. + * \param[in] CoordCorners - Coordiantes of the corners. + * \param[in] DShapeFunction - Shape function information + */ + su2double ShapeFunc_Quadrilateral(su2double Xi, su2double Eta, su2double CoordCorners[8][3], + su2double DShapeFunction[8][4]); + + /*! + * \brief Check the domain points vertex that are going to be moved. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetDomainDisplacements(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 derivatives of the boundary nodes. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetBoundaryDerivatives(CGeometry* geometry, CConfig* config, bool ForwardProjectionDerivative); + + /*! + * \brief Store the number of iterations when moving the mesh. + * \param[in] val_nIterMesh - Number of iterations. + */ + inline void Set_nIterMesh(unsigned long val_nIterMesh) { nIterMesh = val_nIterMesh; } + + /*! + * \brief Retrieve the number of iterations when moving the mesh. + * \param[out] Number of iterations. + */ + inline unsigned long Get_nIterMesh() const { return nIterMesh; } +}; + + diff --git a/Common/include/grid_movement/CRadialBasisFunctionInterpolation.hpp b/Common/include/grid_movement/CRadialBasisFunctionInterpolation.hpp new file mode 100644 index 000000000000..daa99b5059d6 --- /dev/null +++ b/Common/include/grid_movement/CRadialBasisFunctionInterpolation.hpp @@ -0,0 +1,247 @@ +/*! + * \file CRadialBasisFunctionInterpolation.hpp + * \brief Headers of the CRadialBasisFunctionInterpolation class. + * \author F. van Steen + * \version 8.0.1 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once +#include "CVolumetricMovement.hpp" +#include "CRadialBasisFunctionNode.hpp" +#include "../../include/toolboxes/CSymmetricMatrix.hpp" + +/*! + * \class CRadialBasisFunctionInterpolation + * \brief Class for moving the volumetric numerical grid using Radial Basis Function interpolation. + * \author F. van Steen + */ + +class CRadialBasisFunctionInterpolation : public CVolumetricMovement { +protected: + + vector* ControlNodes = nullptr; /*!< \brief Vector with control nodes*/ + vector BoundNodes; /*!< \brief Vector with boundary nodes.*/ + vector ReducedControlNodes; /*!< \brief Vector with selected control nodes in data reduction algorithm. */ + + + vector CtrlNodeDeformation; /*!< \brief Control Node Deformation.*/ + vector InterpCoeff; /*!< \brief Control node interpolation coefficients.*/ + + unsigned long nCtrlNodesGlobal{0}; /*!< \brief Total number of control nodes.*/ + su2activematrix CtrlCoords; /*!< \brief Coordinates of the control nodes.*/ + + su2double MaxErrorGlobal{0.0}; /*!< \brief Maximum error data reduction algorithm.*/ + + +public: + + /*! + * \brief Constructor of the class. + */ + CRadialBasisFunctionInterpolation(CGeometry* geometry, CConfig* config); + + /*! + * \brief Destructor of the class. + */ + ~CRadialBasisFunctionInterpolation(void) override; + + /*! + * \brief Grid deformation using the spring analogy method. + * \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(CGeometry* geometry, CConfig* config, bool UpdateGeo, bool Derivative, + bool ForwardProjectionDerivative); + + /*! + * \brief Selecting unique set of boundary nodes based on marker information. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetBoundNodes(CGeometry* geometry, CConfig* config); + + /*! + * \brief Selecting internal nodes for the volumetric deformation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] internalNode - Internal nodes. + */ + void SetInternalNodes(CGeometry* geometry, CConfig* config, vector& internalNodes); + + /*! + * \brief Assigning the control nodes. + * \param[in] config -Definition of the particular problem. + * */ + + void SetCtrlNodes(CConfig* config); + + /*! + * \brief Solving the RBF system to obtain the interpolation coefficients. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] type - Type of radial basis function. + * \param[in] radius - Support radius of the radial basis function. + */ + + void SolveRBF_System(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius); + + /*! + * \brief Obtaining the interpolation coefficients of the control nodes. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] type - Type of radial basis function. + * \param[in] radius - Support radius of the radial basis function. + */ + + void GetInterpCoeffs(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius); + + + /*! + * \brief Gathering of all control node coordinates. + * \param[in] geometry - Geometrical definition of the problem. + */ + void SetCtrlNodeCoords(CGeometry* geometry); + + /*! + * \brief Build the deformation vector with surface displacements of the control nodes. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetDeformation(CGeometry* geometry, CConfig* config); + + /*! + * \brief Computation of the interpolation matrix and inverting in. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] type - Type of radial basis function. + * \param[in] radius - Support radius of the radial basis function. + * \param[in] invInterpMat - Inverse of the interpolation matrix. + */ + void ComputeInterpolationMatrix(CGeometry* geometry, const RADIAL_BASIS& type, const su2double radius, su2passivematrix& invInterpMat); + + /*! + * \brief Computation of interpolation coefficients + * \param[in] invInterpMat - Inverse of interpolation matrix + */ + void ComputeInterpCoeffs(su2passivematrix& invInterpMat); + + /*! + * \brief Finding initial data reduction control nodes based on maximum deformation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] maxErrorNodeLocal - Local maximum error node. + * \param[in] maxErrorLocal - Local maximum error. + */ + void GetInitMaxErrorNode(CGeometry* geometry, CConfig* config, unsigned long& maxErrorNodeLocal, su2double& maxErrorLocal); + + /*! + * \brief Addition of control node to the reduced set. + * \param[in] maxErrorNode - Node with maximum error to be added. + */ + void AddControlNode(unsigned long maxErrorNode); + + /*! + * \brief Compute global number of control nodes. + */ + void Get_nCtrlNodesGlobal(); + + /*! + * \brief Compute interpolation error. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] type - Type of radial basis function. + * \param[in] radius - Support radius of the radial basis function. + * \param[in] maxErrorNodeLocal - Local maximum error node. + * \param[in] maxErrorLocal - Local maximum error. + */ + void GetInterpError(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius, unsigned long& maxErrorNodeLocal, su2double& maxErrorLocal); + + /*! + * \brief Compute error of single node. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] type - Type of radial basis function. + * \param[in] radius - Support radius of the radial basis function. + * \param[in] iNode - Local node in consideration. + * \param[in] localError - Local error. + */ + void GetNodalError(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius, unsigned long iNode, su2double* localError); + + /*! + * \brief Updating the grid coordinates. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] type - Type of radial basis function. + * \param[in] radius - Support radius of the radial basis function. + * \param[in] internalNodes - Internal nodes. + */ + void UpdateGridCoord(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius, const vector& internalNodes); + + /*! + * \brief Updating the internal node coordinates. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] type - Type of radial basis function. + * \param[in] radius - Support radius of the radial basis function. + * \param[in] internalNodes - Internal nodes. + */ + void UpdateInternalCoords(CGeometry* geometry, const RADIAL_BASIS& type, const su2double radius, const vector& internalNodes); + + /*! + * \brief Updating the internal node coordinates. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] type - Type of radial basis function. + * \param[in] radius - Support radius of the radial basis function. + */ + void UpdateBoundCoords(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius); + + /*! + * \brief Apply correction to the nonzero error boundary nodes. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] type - Type of radial basis function. + * \param[in] internalNodes - Internal nodes. + */ + void SetCorrection(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const vector& internalNodes); + + /*! + * \brief Custom comparison function, for sorting the CRadialBasisFunctionNode objects based on their index. + * \param[in] a - First considered Radial Basis Function Node. + * \param[in] b - Second considered Radial Basis Function Node. + * \return True if index of a is smaller than index of b. + */ + inline static bool HasSmallerIndex(CRadialBasisFunctionNode* a, CRadialBasisFunctionNode* b){ + return a->GetIndex() < b->GetIndex(); + } + + /*! + * \brief Custom equality function, for obtaining a unique set of CRadialBasisFunctionNode objects. + * \param[in] a - First considered Radial Basis Function Node. + * \param[in] b - Second considered Radial Basis Function Node. + * \return True if index of a and b are equal. + */ + inline static bool HasEqualIndex(CRadialBasisFunctionNode* a, CRadialBasisFunctionNode* b){ + return a->GetIndex() == b->GetIndex(); + } +}; \ No newline at end of file diff --git a/Common/include/grid_movement/CRadialBasisFunctionNode.hpp b/Common/include/grid_movement/CRadialBasisFunctionNode.hpp new file mode 100644 index 000000000000..0b080b87ad14 --- /dev/null +++ b/Common/include/grid_movement/CRadialBasisFunctionNode.hpp @@ -0,0 +1,89 @@ +/*! + * \file CRadialBasisFunctionNode.hpp + * \brief Declaration of the RBF node class that stores nodal information for the RBF interpolation. + * \author F. van Steen + * \version 8.0.1 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "../../../Common/include/geometry/CGeometry.hpp" + +/*! + * \class CRadialBasisFunctionNode + * \brief Class for defining a Radial Basis Function Node + * \author F. van Steen + */ + +class CRadialBasisFunctionNode{ + protected: + unsigned long idx; /*!< \brief Global index. */ + unsigned short marker_idx; /*!< \brief Marker index. */ + unsigned long vertex_idx; /*!< \brief Vertex index. */ + + su2double error[3]; /*!< \brief Nodal data reduction error; */ + + public: + + /*! + * \brief Constructor of the class. + * \param[in] idx_val - Local node index. + * \param[in] marker_val - Local marker index. + * \param[in] vertex_val - Local vertex index. + */ + CRadialBasisFunctionNode(unsigned long idx_val, unsigned short marker_val, unsigned long vertex_val); + + /*! + * \brief Returns local global index. + * \return Local node index. + */ + inline unsigned long GetIndex(){return idx;} + + /*! + * \brief Returns local vertex index. + * \return Local vertex index. + */ + inline unsigned long GetVertex(){return vertex_idx;} + + /*! + * \brief Returns local marker index. + * \return Local marker index. + */ + inline unsigned short GetMarker(){return marker_idx;} + + /*! + * \brief Set error of the RBF node. + * \param val_error - Nodal error. + * \param nDim - Number of dimensions. + */ + + inline void SetError(const su2double* val_error, unsigned short nDim) { + for (auto iDim = 0u; iDim < nDim; iDim++) error[iDim] = val_error[iDim]; + } + + /*! + * \brief Get nodal error. + * \return Nodal error. + */ + inline su2double* GetError(){ return error;} +}; \ No newline at end of file diff --git a/Common/include/grid_movement/CVolumetricMovement.hpp b/Common/include/grid_movement/CVolumetricMovement.hpp index 95f324bb9024..dcb787894fe7 100644 --- a/Common/include/grid_movement/CVolumetricMovement.hpp +++ b/Common/include/grid_movement/CVolumetricMovement.hpp @@ -28,9 +28,7 @@ #pragma once #include "CGridMovement.hpp" -#include "../linear_algebra/CSysMatrix.hpp" -#include "../linear_algebra/CSysVector.hpp" -#include "../linear_algebra/CSysSolve.hpp" + /*! * \class CVolumetricMovement @@ -39,23 +37,7 @@ */ class CVolumetricMovement : public CGridMovement { 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. +*/ - -#ifndef CODI_FORWARD_TYPE - CSysMatrix StiffMatrix; /*!< \brief Stiffness matrix of the elasticity problem. */ - CSysSolve System; /*!< \brief Linear solver/smoother. */ -#else - CSysMatrix StiffMatrix; - CSysSolve System; -#endif - CSysVector LinSysSol; - CSysVector LinSysRes; + unsigned short nDim; /*!< \brief Number of dimensions. */ public: /*! @@ -64,22 +46,16 @@ class CVolumetricMovement : public CGridMovement { CVolumetricMovement(void); /*! - * \brief Constructor of the class. - */ - CVolumetricMovement(CGeometry* geometry, CConfig* config); + * \brief Constructor of the Class. + */ + + CVolumetricMovement(CGeometry* geometry); /*! * \brief Destructor of the class. */ ~CVolumetricMovement(void) override; - /*! - * \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. @@ -94,106 +70,6 @@ class CVolumetricMovement : public CGridMovement { */ void UpdateMultiGrid(CGeometry** geometry, CConfig* config); - /*! - * \brief Compute the stiffness matrix for grid deformation using spring analogy. - * \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. - */ - su2double SetFEAMethodContributions_Elem(CGeometry* geometry, CConfig* config); - - /*! - * \brief Build the stiffness matrix for a 3-D hexahedron element. The result will be placed in StiffMatrix_Elem. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] StiffMatrix_Elem - Element stiffness matrix to be filled. - * \param[in] CoordCorners - Index value for Node 1 of the current hexahedron. - * \param[in] PointCorners - Index values for element corners - * \param[in] nNodes - Number of nodes defining the element. - * \param[in] scale - */ - void SetFEA_StiffMatrix3D(CGeometry* geometry, CConfig* config, su2double** StiffMatrix_Elem, - unsigned long PointCorners[8], su2double CoordCorners[8][3], unsigned short nNodes, - su2double ElemVolume, su2double ElemDistance); - - /*! - * \brief Build the stiffness matrix for a 3-D hexahedron element. The result will be placed in StiffMatrix_Elem. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] StiffMatrix_Elem - Element stiffness matrix to be filled. - * \param[in] CoordCorners - Index value for Node 1 of the current hexahedron. - * \param[in] PointCorners - Index values for element corners - * \param[in] nNodes - Number of nodes defining the element. - * \param[in] scale - */ - void SetFEA_StiffMatrix2D(CGeometry* geometry, CConfig* config, su2double** StiffMatrix_Elem, - unsigned long PointCorners[8], su2double CoordCorners[8][3], unsigned short nNodes, - su2double ElemVolume, su2double ElemDistance); - - /*! - * \brief Shape functions and derivative of the shape functions - * \param[in] Xi - Local coordinates. - * \param[in] Eta - Local coordinates. - * \param[in] Zeta - Local coordinates. - * \param[in] CoordCorners - Coordiantes of the corners. - * \param[in] DShapeFunction - Shape function information - */ - su2double ShapeFunc_Hexa(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3], - su2double DShapeFunction[8][4]); - - /*! - * \brief Shape functions and derivative of the shape functions - * \param[in] Xi - Local coordinates. - * \param[in] Eta - Local coordinates. - * \param[in] Zeta - Local coordinates. - * \param[in] CoordCorners - Coordiantes of the corners. - * \param[in] DShapeFunction - Shape function information - */ - su2double ShapeFunc_Tetra(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3], - su2double DShapeFunction[8][4]); - - /*! - * \brief Shape functions and derivative of the shape functions - * \param[in] Xi - Local coordinates. - * \param[in] Eta - Local coordinates. - * \param[in] Zeta - Local coordinates. - * \param[in] CoordCorners - Coordiantes of the corners. - * \param[in] DShapeFunction - Shape function information - */ - su2double ShapeFunc_Pyram(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3], - su2double DShapeFunction[8][4]); - - /*! - * \brief Shape functions and derivative of the shape functions - * \param[in] Xi - Local coordinates. - * \param[in] Eta - Local coordinates. - * \param[in] Zeta - Local coordinates. - * \param[in] CoordCorners - Coordiantes of the corners. - * \param[in] DShapeFunction - Shape function information - */ - su2double ShapeFunc_Prism(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3], - su2double DShapeFunction[8][4]); - - /*! - * \brief Shape functions and derivative of the shape functions - * \param[in] Xi - Local coordinates. - * \param[in] Eta - Local coordinates. - * \param[in] CoordCorners - Coordiantes of the corners. - * \param[in] DShapeFunction - Shape function information - */ - su2double ShapeFunc_Triangle(su2double Xi, su2double Eta, su2double CoordCorners[8][3], - su2double DShapeFunction[8][4]); - - /*! - * \brief Shape functions and derivative of the shape functions - * \param[in] Xi - Local coordinates. - * \param[in] Eta - Local coordinates. - * \param[in] CoordCorners - Coordiantes of the corners. - * \param[in] DShapeFunction - Shape function information - */ - su2double ShapeFunc_Quadrilateral(su2double Xi, su2double Eta, su2double CoordCorners[8][3], - su2double DShapeFunction[8][4]); - /*! * \brief Compute the shape functions for hexahedron * \param[in] CoordCorners - coordinates of the cornes of the hexahedron. @@ -230,15 +106,6 @@ class CVolumetricMovement : public CGridMovement { */ su2double GetQuadrilateral_Area(su2double CoordCorners[8][3]) const; - /*! - * \brief Add the stiffness matrix for a 2-D triangular element to the global stiffness matrix for the entire mesh - * (node-based). \param[in] geometry - Geometrical definition of the problem. \param[in] StiffMatrix_Elem - Element - * stiffness matrix to be filled. \param[in] PointCorners - Index values for element corners \param[in] nNodes - - * Number of nodes defining the element. - */ - void AddFEA_StiffMatrix(CGeometry* geometry, su2double** StiffMatrix_Elem, unsigned long PointCorners[8], - unsigned short nNodes); - /*! * \brief Check for negative volumes (all elements) after performing grid deformation. * \param[in] geometry - Geometrical definition of the problem. @@ -253,29 +120,7 @@ class CVolumetricMovement : public CGridMovement { * \param[in] Screen_Output - determines if text is written to screen */ void ComputenNonconvexElements(CGeometry* geometry, bool Screen_Output); - - /*! - * \brief Compute the minimum distance to the nearest solid surface. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void ComputeSolid_Wall_Distance(CGeometry* geometry, CConfig* config, su2double& MinDistance, - su2double& MaxDistance) const; - - /*! - * \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 Check the domain points vertex that are going to be moved. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void SetDomainDisplacements(CGeometry* geometry, CConfig* config); - + /*! * \brief Unsteady grid movement using rigid mesh rotation. * \param[in] geometry - Geometrical definition of the problem. @@ -343,8 +188,8 @@ class CVolumetricMovement : public CGridMovement { * \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(CGeometry* geometry, CConfig* config, bool UpdateGeo, bool Derivative = false, - bool ForwardProjectionDerivative = false); + inline virtual void SetVolume_Deformation(CGeometry* geometry, CConfig* config, bool UpdateGeo, bool Derivative = false, + bool ForwardProjectionDerivative = false){}; /*! * \brief Grid deformation using the spring analogy method. @@ -356,32 +201,6 @@ class CVolumetricMovement : public CGridMovement { inline virtual void SetVolume_Deformation_Elas(CGeometry* geometry, CConfig* config, bool UpdateGeo, bool screen_output, bool Derivative = false) {} - /*! - * \brief Set the derivatives of the boundary nodes. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void SetBoundaryDerivatives(CGeometry* geometry, CConfig* config, bool ForwardProjectionDerivative); - - /*! - * \brief Update the derivatives 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_Derivatives(CGeometry* geometry, CConfig* config, bool ForwardProjectionDerivative); - - /*! - * \brief Store the number of iterations when moving the mesh. - * \param[in] val_nIterMesh - Number of iterations. - */ - inline void Set_nIterMesh(unsigned long val_nIterMesh) { nIterMesh = val_nIterMesh; } - - /*! - * \brief Retrieve the number of iterations when moving the mesh. - * \param[out] Number of iterations. - */ - inline unsigned long Get_nIterMesh() const { return nIterMesh; } - /*! * \brief Set the boundary dependencies in the mesh side of the problem * \param[in] geometry - Geometrical definition of the problem. diff --git a/Common/include/grid_movement/CVolumetricMovementFactory.hpp b/Common/include/grid_movement/CVolumetricMovementFactory.hpp new file mode 100644 index 000000000000..5ce673f1d58a --- /dev/null +++ b/Common/include/grid_movement/CVolumetricMovementFactory.hpp @@ -0,0 +1,42 @@ +/*! + * \file CVolumetricMovementFactory.hpp + * \brief Factory to generate volumetric mover objects. + * \version 8.0.1 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +class CConfig; +class CGeometry; +class CVolumetricMovement; + +namespace CVolumetricMovementFactory { +/*! + * \brief Factory method for CVolumetricMovement objects. + * \param[in] geometry_container - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \return Pointer to the allocated volumetric mover. + */ + + CVolumetricMovement* CreateCVolumetricMovement(CGeometry* geometry, CConfig* config); +} \ No newline at end of file diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 7199a310193d..662da7b39389 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2642,6 +2642,19 @@ static const MapType Sobolev_Modus_Map = { MakePair("ONLY_GRADIENT", ENUM_SOBOLEV_MODUS::ONLY_GRAD) }; +/*! + * \brief Type of mesh deformation + */ +enum class DEFORM_KIND { + ELASTIC, /*!< \brief Linear elasticity method. */ + RBF /*!< \brief Radial basis function interpolation. */ +}; +static const MapType Deform_Kind_Map = { + MakePair("ELASTIC", DEFORM_KIND::ELASTIC) + MakePair("RBF", DEFORM_KIND::RBF) +}; + + #undef MakePair /* END_CONFIG_ENUMS */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 18da76556c40..742cec397f03 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -837,6 +837,7 @@ void CConfig::SetPointersNull() { Marker_CfgFile_ZoneInterface = nullptr; Marker_CfgFile_Deform_Mesh = nullptr; Marker_All_Deform_Mesh = nullptr; Marker_CfgFile_Deform_Mesh_Sym_Plane = nullptr; Marker_All_Deform_Mesh_Sym_Plane = nullptr; + Marker_CfgFile_Deform_Mesh_Internal = nullptr; Marker_All_Deform_Mesh_Internal = nullptr; Marker_CfgFile_Fluid_Load = nullptr; Marker_All_Fluid_Load = nullptr; Marker_CfgFile_SobolevBC = nullptr; Marker_All_SobolevBC = nullptr; @@ -863,7 +864,8 @@ void CConfig::SetPointersNull() { Marker_Euler = nullptr; Marker_FarField = nullptr; Marker_Custom = nullptr; Marker_SymWall = nullptr; Marker_PerBound = nullptr; Marker_PerDonor = nullptr; Marker_NearFieldBound = nullptr; Marker_Inlet_Turb = nullptr; - Marker_Deform_Mesh = nullptr; Marker_Deform_Mesh_Sym_Plane= nullptr; Marker_Fluid_Load = nullptr; + Marker_Deform_Mesh = nullptr; Marker_Deform_Mesh_Sym_Plane= nullptr; Marker_Deform_Mesh_Internal= nullptr; + Marker_Fluid_Load = nullptr; Marker_Inlet = nullptr; Marker_Outlet = nullptr; Marker_Inlet_Species = nullptr; Marker_Supersonic_Inlet = nullptr; Marker_Supersonic_Outlet = nullptr; Marker_Smoluchowski_Maxwell= nullptr; Marker_Isothermal = nullptr; Marker_HeatFlux = nullptr; Marker_EngineInflow = nullptr; @@ -1497,6 +1499,8 @@ void CConfig::SetConfig_Options() { addStringListOption("MARKER_DEFORM_MESH", nMarker_Deform_Mesh, Marker_Deform_Mesh); /*!\brief MARKER_DEFORM_MESH_SYM_PLANE\n DESCRIPTION: Symmetry plane for mesh deformation only \ingroup Config*/ addStringListOption("MARKER_DEFORM_MESH_SYM_PLANE", nMarker_Deform_Mesh_Sym_Plane, Marker_Deform_Mesh_Sym_Plane); + /*!\brief MARKER_DEFORM_MESH_INTERNAL\n DESCRIPTION: Internal marker for mesh deformation only \ingroup Config*/ + addStringListOption("MARKER_DEFORM_MESH_INTERNAL", nMarker_Deform_Mesh_Internal, Marker_Deform_Mesh_Internal); /*!\brief MARKER_FLUID_LOAD\n DESCRIPTION: Marker(s) in which the flow load is computed/applied \ingroup Config*/ addStringListOption("MARKER_FLUID_LOAD", nMarker_Fluid_Load, Marker_Fluid_Load); /*!\brief MARKER_FSI_INTERFACE \n DESCRIPTION: ZONE interface boundary marker(s) \ingroup Config*/ @@ -2375,6 +2379,14 @@ void CConfig::SetConfig_Options() { addDoubleOption("DEFORM_LINEAR_SOLVER_ERROR", Deform_Linear_Solver_Error, 1E-14); /* DESCRIPTION: Maximum number of iterations of the linear solver for the implicit formulation */ addUnsignedLongOption("DEFORM_LINEAR_SOLVER_ITER", Deform_Linear_Solver_Iter, 1000); + /* DESCRIPTION: Type of mesh deformation */ + addEnumOption("DEFORM_KIND", Deform_Kind, Deform_Kind_Map, DEFORM_KIND::ELASTIC); + /* DESCRIPTION: Use of data reduction methods for RBF interpolated mesh deformation */ + addBoolOption("RBF_DATA_REDUCTION", RBF_DataReduction, true); + /* DESCRIPTION: Tolerance for the data reduction methods used in RBF mesh deformation. */ + addDoubleOption("RBF_GREEDY_TOLERANCE", RBF_GreedyTolerance, 1E-2); + /* DESCRIPTION: Tolerance for the data reduction methods used in RBF mesh deformation. */ + addDoubleOption("RBF_GREEDY_CORRECTION_FACTOR", RBF_GreedyCorrectionFactor, 1E-2); /*!\par CONFIG_CATEGORY: FEM flow solver definition \ingroup Config*/ /*--- Options related to the finite element flow solver---*/ @@ -5563,7 +5575,7 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { iMarker_Monitoring, iMarker_Designing, iMarker_GeoEval, iMarker_Plotting, iMarker_Analyze, iMarker_DV, iMarker_Moving, iMarker_SobolevBC, iMarker_PyCustom, iMarker_Supersonic_Inlet, iMarker_Supersonic_Outlet, iMarker_Clamped, iMarker_ZoneInterface, iMarker_CHTInterface, iMarker_Load_Dir, iMarker_Disp_Dir, - iMarker_Fluid_Load, iMarker_Deform_Mesh, iMarker_Deform_Mesh_Sym_Plane, + iMarker_Fluid_Load, iMarker_Deform_Mesh, iMarker_Deform_Mesh_Sym_Plane, iMarker_Deform_Mesh_Internal, iMarker_ActDiskInlet, iMarker_ActDiskOutlet, iMarker_Turbomachinery, iMarker_MixingPlaneInterface; @@ -5607,6 +5619,7 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { Marker_All_Moving = new unsigned short[nMarker_All] (); // Store whether the boundary should be in motion. Marker_All_Deform_Mesh = new unsigned short[nMarker_All] (); // Store whether the boundary is deformable. Marker_All_Deform_Mesh_Sym_Plane = new unsigned short[nMarker_All] (); //Store wheter the boundary will follow the deformation + Marker_All_Deform_Mesh_Internal = new unsigned short[nMarker_All] (); //Store wheter the boundary will follow the deformation Marker_All_Fluid_Load = new unsigned short[nMarker_All] (); // Store whether the boundary computes/applies fluid loads. Marker_All_PyCustom = new unsigned short[nMarker_All] (); // Store whether the boundary is Python customizable. Marker_All_PerBound = new short[nMarker_All] (); // Store whether the boundary belongs to a periodic boundary. @@ -5633,6 +5646,7 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { Marker_CfgFile_Moving = new unsigned short[nMarker_CfgFile] (); Marker_CfgFile_Deform_Mesh = new unsigned short[nMarker_CfgFile] (); Marker_CfgFile_Deform_Mesh_Sym_Plane= new unsigned short[nMarker_CfgFile] (); + Marker_CfgFile_Deform_Mesh_Internal = new unsigned short[nMarker_CfgFile] (); Marker_CfgFile_Fluid_Load = new unsigned short[nMarker_CfgFile] (); Marker_CfgFile_PerBound = new unsigned short[nMarker_CfgFile] (); Marker_CfgFile_Turbomachinery = new unsigned short[nMarker_CfgFile] (); @@ -6047,6 +6061,13 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) { Marker_CfgFile_Deform_Mesh_Sym_Plane[iMarker_CfgFile] = YES; } + for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) { + Marker_CfgFile_Deform_Mesh_Internal[iMarker_CfgFile] = NO; + for (iMarker_Deform_Mesh_Internal = 0; iMarker_Deform_Mesh_Internal < nMarker_Deform_Mesh_Internal; iMarker_Deform_Mesh_Internal++) + if (Marker_CfgFile_TagBound[iMarker_CfgFile] == Marker_Deform_Mesh_Internal[iMarker_Deform_Mesh_Internal]) + Marker_CfgFile_Deform_Mesh_Internal[iMarker_CfgFile] = YES; + } + for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) { Marker_CfgFile_Fluid_Load[iMarker_CfgFile] = NO; for (iMarker_Fluid_Load = 0; iMarker_Fluid_Load < nMarker_Fluid_Load; iMarker_Fluid_Load++) @@ -6075,8 +6096,8 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { unsigned short iMarker_Euler, iMarker_Custom, iMarker_FarField, iMarker_SymWall, iMarker_PerBound, iMarker_NearFieldBound, iMarker_Fluid_InterfaceBound, iMarker_Inlet, iMarker_Riemann, - iMarker_Deform_Mesh, iMarker_Deform_Mesh_Sym_Plane, iMarker_Fluid_Load, - iMarker_Smoluchowski_Maxwell, iWall_Catalytic, + iMarker_Deform_Mesh, iMarker_Deform_Mesh_Sym_Plane, iMarker_Deform_Mesh_Internal, + iMarker_Fluid_Load, iMarker_Smoluchowski_Maxwell, iWall_Catalytic, iMarker_Giles, iMarker_Outlet, iMarker_Isothermal, iMarker_HeatFlux, iMarker_HeatTransfer, iMarker_EngineInflow, iMarker_EngineExhaust, iMarker_Displacement, iMarker_Damper, iMarker_Load, iMarker_Internal, iMarker_Monitoring, @@ -6457,7 +6478,19 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { if (val_software == SU2_COMPONENT::SU2_DEF) { cout << endl <<"---------------- Grid deformation parameters ( Zone " << iZone << " ) ----------------" << endl; - cout << "Grid deformation using a linear elasticity method." << endl; + cout << "Grid deformation using "; + if (Deform_Kind == DEFORM_KIND::RBF){ + cout << "Radial Basis Function interpolation method.\nRadial Basis Function: "; + switch(Kind_RadialBasisFunction){ + case RADIAL_BASIS::WENDLAND_C2: cout << "Wendland C2." << endl; break; + case RADIAL_BASIS::INV_MULTI_QUADRIC: cout << "inversed multi quartic biharmonic spline." << endl; break; + case RADIAL_BASIS::GAUSSIAN: cout << "Guassian." << endl; break; + case RADIAL_BASIS::THIN_PLATE_SPLINE: cout << "thin plate spline." << endl; break; + case RADIAL_BASIS::MULTI_QUADRIC: cout << "multi quartic biharmonic spline." << endl; break; + } + }else{ + cout << "linear elasticity method." << endl; + } if (Hold_GridFixed == YES) cout << "Hold some regions of the mesh fixed (hardcode implementation)." << endl; } @@ -7349,6 +7382,15 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { BoundaryTable.PrintFooter(); } + if (nMarker_Deform_Mesh_Internal!= 0) { + BoundaryTable << "Internal freely deformable mesh boundary"; + for (iMarker_Deform_Mesh_Internal = 0; iMarker_Deform_Mesh_Internal < nMarker_Deform_Mesh_Internal; iMarker_Deform_Mesh_Internal++) { + BoundaryTable << Marker_Deform_Mesh_Internal[iMarker_Deform_Mesh_Internal]; + if (iMarker_Deform_Mesh_Internal < nMarker_Deform_Mesh_Internal-1) BoundaryTable << " "; + } + BoundaryTable.PrintFooter(); + } + if (nMarker_Fluid_Load != 0) { BoundaryTable << "Fluid loads boundary"; for (iMarker_Fluid_Load = 0; iMarker_Fluid_Load < nMarker_Fluid_Load; iMarker_Fluid_Load++) { @@ -7903,6 +7945,13 @@ unsigned short CConfig::GetMarker_CfgFile_Deform_Mesh_Sym_Plane(const string& va return Marker_CfgFile_Deform_Mesh_Sym_Plane[iMarker_CfgFile]; } +unsigned short CConfig::GetMarker_CfgFile_Deform_Mesh_Internal(const string& val_marker) const { + unsigned short iMarker_CfgFile; + for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) + if (Marker_CfgFile_TagBound[iMarker_CfgFile] == val_marker) break; + return Marker_CfgFile_Deform_Mesh_Internal[iMarker_CfgFile]; +} + unsigned short CConfig::GetMarker_CfgFile_Fluid_Load(const string& val_marker) const { unsigned short iMarker_CfgFile; for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) @@ -8047,6 +8096,9 @@ CConfig::~CConfig() { delete[] Marker_CfgFile_Deform_Mesh_Sym_Plane; delete[] Marker_All_Deform_Mesh_Sym_Plane; + delete[] Marker_CfgFile_Deform_Mesh_Internal; + delete[] Marker_All_Deform_Mesh_Internal; + delete[] Marker_CfgFile_Fluid_Load; delete[] Marker_All_Fluid_Load; @@ -8839,6 +8891,15 @@ unsigned short CConfig::GetMarker_Deform_Mesh_Sym_Plane(const string& val_marker return iMarker; } +unsigned short CConfig::GetMarker_Deform_Mesh_Internal(const string& val_marker) const { + unsigned short iMarker; + + /*--- Find the marker for this internal boundary. ---*/ + for (iMarker = 0; iMarker < nMarker_Deform_Mesh_Internal; iMarker++) + if (Marker_Deform_Mesh_Internal[iMarker] == val_marker) break; + return iMarker; +} + unsigned short CConfig::GetMarker_Fluid_Load(const string& val_marker) const { unsigned short iMarker_Fluid_Load; diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 8d245ce1206c..c64499b817ef 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -3365,6 +3365,7 @@ void CPhysicalGeometry::SetBoundaries(CConfig* config) { config->SetMarker_All_Moving(iMarker, config->GetMarker_CfgFile_Moving(Marker_Tag)); config->SetMarker_All_Deform_Mesh(iMarker, config->GetMarker_CfgFile_Deform_Mesh(Marker_Tag)); config->SetMarker_All_Deform_Mesh_Sym_Plane(iMarker, config->GetMarker_CfgFile_Deform_Mesh_Sym_Plane(Marker_Tag)); + config->SetMarker_All_Deform_Mesh_Internal(iMarker, config->GetMarker_CfgFile_Deform_Mesh_Internal(Marker_Tag)); config->SetMarker_All_Fluid_Load(iMarker, config->GetMarker_CfgFile_Fluid_Load(Marker_Tag)); config->SetMarker_All_PyCustom(iMarker, config->GetMarker_CfgFile_PyCustom(Marker_Tag)); config->SetMarker_All_PerBound(iMarker, config->GetMarker_CfgFile_PerBound(Marker_Tag)); @@ -3389,6 +3390,7 @@ void CPhysicalGeometry::SetBoundaries(CConfig* config) { config->SetMarker_All_Moving(iMarker, NO); config->SetMarker_All_Deform_Mesh(iMarker, NO); config->SetMarker_All_Deform_Mesh_Sym_Plane(iMarker, NO); + config->SetMarker_All_Deform_Mesh_Internal(iMarker, NO); config->SetMarker_All_Fluid_Load(iMarker, NO); config->SetMarker_All_PyCustom(iMarker, NO); config->SetMarker_All_PerBound(iMarker, NO); @@ -3758,6 +3760,7 @@ void CPhysicalGeometry::LoadUnpartitionedSurfaceElements(CConfig* config, CMeshR config->SetMarker_All_Moving(iMarker, config->GetMarker_CfgFile_Moving(Marker_Tag)); config->SetMarker_All_Deform_Mesh(iMarker, config->GetMarker_CfgFile_Deform_Mesh(Marker_Tag)); config->SetMarker_All_Deform_Mesh_Sym_Plane(iMarker, config->GetMarker_CfgFile_Deform_Mesh_Sym_Plane(Marker_Tag)); + config->SetMarker_All_Deform_Mesh_Internal(iMarker, config->GetMarker_CfgFile_Deform_Mesh_Internal(Marker_Tag)); config->SetMarker_All_Fluid_Load(iMarker, config->GetMarker_CfgFile_Fluid_Load(Marker_Tag)); config->SetMarker_All_PyCustom(iMarker, config->GetMarker_CfgFile_PyCustom(Marker_Tag)); config->SetMarker_All_PerBound(iMarker, config->GetMarker_CfgFile_PerBound(Marker_Tag)); diff --git a/Common/src/grid_movement/CLinearElasticity.cpp b/Common/src/grid_movement/CLinearElasticity.cpp new file mode 100644 index 000000000000..1f7b6f99a55b --- /dev/null +++ b/Common/src/grid_movement/CLinearElasticity.cpp @@ -0,0 +1,1479 @@ +/*! + * \file CLinearElasticity.cpp + * \brief Subroutines for moving mesh volume elements using the linear elasticity analogy + * \author F. Palacios, T. Economon, S. Padron + * \version 8.0.1 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/grid_movement/CLinearElasticity.hpp" +#include "../../include/adt/CADTPointsOnlyClass.hpp" + +CLinearElasticity::CLinearElasticity(CGeometry* geometry, CConfig* config):CVolumetricMovement(geometry),System(LINEAR_SOLVER_MODE::MESH_DEFORM){ + + /*--- Initialize the number of spatial dimensions, length of the state + vector (same as spatial dimensions for grid deformation), and grid nodes. ---*/ + + nVar = geometry->GetnDim(); + nPoint = geometry->GetnPoint(); + nPointDomain = geometry->GetnPointDomain(); + + nIterMesh = 0; + + /*--- Initialize matrix, solution, and r.h.s. structures for the linear solver. ---*/ + if (config->GetVolumetric_Movement() || config->GetSmoothGradient()) { + LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); + LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); + StiffMatrix.Initialize(nPoint, nPointDomain, nVar, nVar, false, geometry, config); + } +}; + +CLinearElasticity::~CLinearElasticity(void) = default; + +void CLinearElasticity::SetVolume_Deformation(CGeometry* geometry, CConfig* config, bool UpdateGeo, bool Derivative, + bool ForwardProjectionDerivative) { + unsigned long Tot_Iter = 0; + su2double MinVolume, MaxVolume; + + /*--- Retrieve number or iterations, tol, output, etc. from config ---*/ + + auto Screen_Output = config->GetDeform_Output(); + auto Nonlinear_Iter = config->GetGridDef_Nonlinear_Iter(); + + /*--- Disable the screen output if we're running SU2_CFD ---*/ + + if (config->GetKind_SU2() == SU2_COMPONENT::SU2_CFD && !Derivative) Screen_Output = false; + if (config->GetSmoothGradient()) Screen_Output = true; + + /*--- Set the number of nonlinear iterations to 1 if Derivative computation is enabled ---*/ + + if (Derivative) Nonlinear_Iter = 1; + + /*--- Loop over the total number of grid deformation iterations. The surface +deformation can be divided into increments to help with stability. In +particular, the linear elasticity equations hold only for small deformations. ---*/ + + for (auto iNonlinear_Iter = 0ul; iNonlinear_Iter < Nonlinear_Iter; iNonlinear_Iter++) { + /*--- Initialize vector and sparse matrix ---*/ + LinSysSol.SetValZero(); + LinSysRes.SetValZero(); + StiffMatrix.SetValZero(); + + /*--- Compute the stiffness matrix entries for all nodes/elements in the + mesh. FEA uses a finite element method discretization of the linear + elasticity equations (transfers element stiffnesses to point-to-point). ---*/ + + MinVolume = SetFEAMethodContributions_Elem(geometry, config); + + /*--- Set the boundary and volume displacements (as prescribed by the + design variable perturbations controlling the surface shape) + as a Dirichlet BC. ---*/ + + SetBoundaryDisplacements(geometry, config); + + /*--- Fix the location of any points in the domain, if requested. ---*/ + + SetDomainDisplacements(geometry, config); + + /*--- Set the boundary derivatives (overrides the actual displacements) ---*/ + + if (Derivative) { + SetBoundaryDerivatives(geometry, config, ForwardProjectionDerivative); + } + + /*--- Communicate any prescribed boundary displacements via MPI, + so that all nodes have the same solution and r.h.s. entries + across all partitions. ---*/ + + CSysMatrixComms::Initiate(LinSysSol, geometry, config); + CSysMatrixComms::Complete(LinSysSol, geometry, config); + + CSysMatrixComms::Initiate(LinSysRes, geometry, config); + CSysMatrixComms::Complete(LinSysRes, geometry, config); + + /*--- Definition of the preconditioner matrix vector multiplication, and linear solver ---*/ + + /*--- To keep legacy behavior ---*/ + System.SetToleranceType(LinearToleranceType::RELATIVE); + + /*--- If we want no derivatives or the direct derivatives, we solve the system using the + * normal matrix vector product and preconditioner. For the mesh sensitivities using + * the discrete adjoint method we solve the system using the transposed matrix. ---*/ + if (!Derivative || ((config->GetKind_SU2() == SU2_COMPONENT::SU2_CFD) && Derivative) || + (config->GetSmoothGradient() && ForwardProjectionDerivative)) { + Tot_Iter = System.Solve(StiffMatrix, LinSysRes, LinSysSol, geometry, config); + } else if (Derivative && (config->GetKind_SU2() == SU2_COMPONENT::SU2_DOT)) { + Tot_Iter = System.Solve_b(StiffMatrix, LinSysRes, LinSysSol, geometry, config); + } + su2double Residual = System.GetResidual(); + + /*--- Update the grid coordinates and cell volumes using the solution + of the linear system (usol contains the x, y, z displacements). ---*/ + + if (!Derivative) { + UpdateGridCoord(geometry, config); + } else { + UpdateGridCoord_Derivatives(geometry, config, ForwardProjectionDerivative); + } + if (UpdateGeo) { + UpdateDualGrid(geometry, config); + } + + + if (!Derivative) { + /*--- Check for failed deformation (negative volumes). ---*/ + + ComputeDeforming_Element_Volume(geometry, MinVolume, MaxVolume, Screen_Output); + + /*--- Calculate amount of nonconvex elements ---*/ + + ComputenNonconvexElements(geometry, Screen_Output); + } + + /*--- Set number of iterations in the mesh update. ---*/ + + Set_nIterMesh(Tot_Iter); + + if (rank == MASTER_NODE && Screen_Output) { + cout << "Non-linear iter.: " << iNonlinear_Iter + 1 << "/" << Nonlinear_Iter << ". Linear iter.: " << Tot_Iter + << ". "; + if (nDim == 2) + cout << "Min. area: " << MinVolume << ". Error: " << Residual << "." << endl; + else + cout << "Min. volume: " << MinVolume << ". Error: " << Residual << "." << endl; + } + } +} + +void CLinearElasticity::UpdateGridCoord(CGeometry* geometry, CConfig* config) { + cout << "updating the grid coordinates" << endl; + unsigned short iDim; + unsigned long iPoint, total_index; + su2double new_coord; + + /*--- Update the grid coordinates using the solution of the linear system +after grid deformation (LinSysSol contains the x, y, z displacements). ---*/ + + for (iPoint = 0; iPoint < nPoint; iPoint++) + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + new_coord = geometry->nodes->GetCoord(iPoint, iDim) + LinSysSol[total_index]; + if (fabs(new_coord) < EPS * EPS) new_coord = 0.0; + geometry->nodes->SetCoord(iPoint, iDim, new_coord); + } + + /*--- LinSysSol contains the non-transformed displacements in the periodic halo cells. + * Hence we still need a communication of the transformed coordinates, otherwise periodicity + * is not maintained. ---*/ + + geometry->InitiateComms(geometry, config, COORDINATES); + geometry->CompleteComms(geometry, config, COORDINATES); +} + +void CLinearElasticity::UpdateGridCoord_Derivatives(CGeometry* geometry, CConfig* config, + bool ForwardProjectionDerivative) { + unsigned short iDim, iMarker; + unsigned long iPoint, total_index, iVertex; + auto* new_coord = new su2double[3]; + + SU2_COMPONENT Kind_SU2 = config->GetKind_SU2(); + + /*--- Update derivatives of the grid coordinates using the solution of the linear system + after grid deformation (LinSysSol contains the derivatives of the x, y, z displacements). ---*/ + if ((config->GetDirectDiff() == D_DESIGN) && (Kind_SU2 == SU2_COMPONENT::SU2_CFD)) { + for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + new_coord[0] = 0.0; + new_coord[1] = 0.0; + new_coord[2] = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + new_coord[iDim] = geometry->nodes->GetCoord(iPoint, iDim); + SU2_TYPE::SetDerivative(new_coord[iDim], SU2_TYPE::GetValue(LinSysSol[total_index])); + } + geometry->nodes->SetCoord(iPoint, new_coord); + } + } else if ((Kind_SU2 == SU2_COMPONENT::SU2_DOT) && !ForwardProjectionDerivative) { + // need to reset here, since we read out the whole vector, but are only interested in boundary derivatives. + if (config->GetSmoothGradient()) { + for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + geometry->SetSensitivity(iPoint, iDim, 0.0); + } + } + } + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if (config->GetSolid_Wall(iMarker) || (config->GetMarker_All_DV(iMarker) == YES)) { + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + if (geometry->nodes->GetDomain(iPoint)) { + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + geometry->SetSensitivity(iPoint, iDim, LinSysSol[total_index]); + } + } + } + } + } + } else if (config->GetSmoothGradient() && ForwardProjectionDerivative) { + for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + geometry->SetSensitivity(iPoint, iDim, LinSysSol[total_index]); + } + } + } + + delete[] new_coord; +} + + +void CLinearElasticity::ComputeSolid_Wall_Distance(CGeometry* geometry, CConfig* config, su2double& MinDistance, + su2double& MaxDistance) const { + unsigned long nVertex_SolidWall, ii, jj, iVertex, iPoint, pointID; + unsigned short iMarker, iDim; + su2double dist, MaxDistance_Local, MinDistance_Local; + int rankID; + + /*--- Initialize min and max distance ---*/ + + MaxDistance = -1E22; + MinDistance = 1E22; + + /*--- Compute the total number of nodes on no-slip boundaries ---*/ + + nVertex_SolidWall = 0; + for (iMarker = 0; iMarker < config->GetnMarker_All(); ++iMarker) { + if (config->GetSolid_Wall(iMarker)) nVertex_SolidWall += geometry->GetnVertex(iMarker); + } + + /*--- Allocate the vectors to hold boundary node coordinates + and its local ID. ---*/ + + vector Coord_bound(nDim * nVertex_SolidWall); + vector PointIDs(nVertex_SolidWall); + + /*--- Retrieve and store the coordinates of the no-slip boundary nodes + and their local point IDs. ---*/ + + ii = 0; + jj = 0; + for (iMarker = 0; iMarker < config->GetnMarker_All(); ++iMarker) { + if (config->GetSolid_Wall(iMarker)) { + for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); ++iVertex) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + PointIDs[jj++] = iPoint; + for (iDim = 0; iDim < nDim; ++iDim) Coord_bound[ii++] = geometry->nodes->GetCoord(iPoint, iDim); + } + } + } + + /*--- Build the ADT of the boundary nodes. ---*/ + + CADTPointsOnlyClass WallADT(nDim, nVertex_SolidWall, Coord_bound.data(), PointIDs.data(), true); + + /*--- Loop over all interior mesh nodes and compute the distances to each + of the no-slip boundary nodes. Store the minimum distance to the wall + for each interior mesh node. ---*/ + + if (WallADT.IsEmpty()) { + /*--- No solid wall boundary nodes in the entire mesh. + Set the wall distance to zero for all nodes. ---*/ + + for (iPoint = 0; iPoint < geometry->GetnPoint(); ++iPoint) geometry->nodes->SetWall_Distance(iPoint, 0.0); + } else { + /*--- Solid wall boundary nodes are present. Compute the wall + distance for all nodes. ---*/ + + for (iPoint = 0; iPoint < geometry->GetnPoint(); ++iPoint) { + WallADT.DetermineNearestNode(geometry->nodes->GetCoord(iPoint), dist, pointID, rankID); + geometry->nodes->SetWall_Distance(iPoint, dist); + + MaxDistance = max(MaxDistance, dist); + + /*--- To discard points on the surface we use > EPS ---*/ + + if (sqrt(dist) > EPS) MinDistance = min(MinDistance, dist); + } + + MaxDistance_Local = MaxDistance; + MaxDistance = 0.0; + MinDistance_Local = MinDistance; + MinDistance = 0.0; + +#ifdef HAVE_MPI + SU2_MPI::Allreduce(&MaxDistance_Local, &MaxDistance, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&MinDistance_Local, &MinDistance, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm()); +#else + MaxDistance = MaxDistance_Local; + MinDistance = MinDistance_Local; +#endif + } +} + +su2double CLinearElasticity::SetFEAMethodContributions_Elem(CGeometry* geometry, CConfig* config) { + unsigned short iVar, iDim, nNodes = 0, iNodes, StiffMatrix_nElem = 0; + unsigned long iElem, PointCorners[8]; + su2double **StiffMatrix_Elem = nullptr, CoordCorners[8][3]; + su2double MinVolume = 0.0, MaxVolume = 0.0, MinDistance = 0.0, MaxDistance = 0.0, ElemVolume = 0.0, + ElemDistance = 0.0; + + bool Screen_Output = config->GetDeform_Output(); + + /*--- Allocate maximum size (quadrilateral and hexahedron) ---*/ + + if (nDim == 2) + StiffMatrix_nElem = 8; + else + StiffMatrix_nElem = 24; + + StiffMatrix_Elem = new su2double*[StiffMatrix_nElem]; + for (iVar = 0; iVar < StiffMatrix_nElem; iVar++) StiffMatrix_Elem[iVar] = new su2double[StiffMatrix_nElem]; + + /*--- Compute min volume in the entire mesh. ---*/ + + ComputeDeforming_Element_Volume(geometry, MinVolume, MaxVolume, Screen_Output); + if (rank == MASTER_NODE && Screen_Output) + cout << "Min. volume: " << MinVolume << ", max. volume: " << MaxVolume << "." << endl; + + /*--- Compute the distance to the nearest surface if needed + as part of the stiffness calculation.. ---*/ + + if ((config->GetDeform_Stiffness_Type() == SOLID_WALL_DISTANCE) || (config->GetDeform_Limit() < 1E6)) { + ComputeSolid_Wall_Distance(geometry, config, MinDistance, MaxDistance); + if (rank == MASTER_NODE && Screen_Output) + cout << "Min. distance: " << MinDistance << ", max. distance: " << MaxDistance << "." << endl; + } + + /*--- Compute contributions from each element by forming the stiffness matrix (FEA) ---*/ + + for (iElem = 0; iElem < geometry->GetnElem(); iElem++) { + if (geometry->elem[iElem]->GetVTK_Type() == TRIANGLE) nNodes = 3; + if (geometry->elem[iElem]->GetVTK_Type() == QUADRILATERAL) nNodes = 4; + if (geometry->elem[iElem]->GetVTK_Type() == TETRAHEDRON) nNodes = 4; + if (geometry->elem[iElem]->GetVTK_Type() == PYRAMID) nNodes = 5; + if (geometry->elem[iElem]->GetVTK_Type() == PRISM) nNodes = 6; + if (geometry->elem[iElem]->GetVTK_Type() == HEXAHEDRON) nNodes = 8; + + for (iNodes = 0; iNodes < nNodes; iNodes++) { + PointCorners[iNodes] = geometry->elem[iElem]->GetNode(iNodes); + for (iDim = 0; iDim < nDim; iDim++) { + CoordCorners[iNodes][iDim] = geometry->nodes->GetCoord(PointCorners[iNodes], iDim); + } + } + + /*--- Extract Element volume and distance to compute the stiffness ---*/ + + ElemVolume = geometry->elem[iElem]->GetVolume(); + + if ((config->GetDeform_Stiffness_Type() == SOLID_WALL_DISTANCE)) { + ElemDistance = 0.0; + for (iNodes = 0; iNodes < nNodes; iNodes++) + ElemDistance += geometry->nodes->GetWall_Distance(PointCorners[iNodes]); + ElemDistance = ElemDistance / (su2double)nNodes; + } + + if (nDim == 2) + SetFEA_StiffMatrix2D(geometry, config, StiffMatrix_Elem, PointCorners, CoordCorners, nNodes, ElemVolume, + ElemDistance); + if (nDim == 3) + SetFEA_StiffMatrix3D(geometry, config, StiffMatrix_Elem, PointCorners, CoordCorners, nNodes, ElemVolume, + ElemDistance); + + AddFEA_StiffMatrix(geometry, StiffMatrix_Elem, PointCorners, nNodes); + } + + /*--- Deallocate memory and exit ---*/ + + for (iVar = 0; iVar < StiffMatrix_nElem; iVar++) delete[] StiffMatrix_Elem[iVar]; + delete[] StiffMatrix_Elem; + + return MinVolume; +} + + +void CLinearElasticity::SetFEA_StiffMatrix2D(CGeometry* geometry, CConfig* config, su2double** StiffMatrix_Elem, + unsigned long PointCorners[8], su2double CoordCorners[8][3], + unsigned short nNodes, su2double ElemVolume, su2double ElemDistance) { + su2double B_Matrix[3][8], D_Matrix[3][3], Aux_Matrix[8][3]; + su2double Xi = 0.0, Eta = 0.0, Det = 0.0, E = 1 / EPS, Lambda = 0.0, Mu = 0.0, Nu = 0.0; + unsigned short iNode, iVar, jVar, kVar, iGauss, nGauss = 0; + su2double DShapeFunction[8][4] = {{0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}}; + su2double Location[4][3], Weight[4]; + unsigned short nVar = geometry->GetnDim(); + + for (iVar = 0; iVar < nNodes * nVar; iVar++) { + for (jVar = 0; jVar < nNodes * nVar; jVar++) { + StiffMatrix_Elem[iVar][jVar] = 0.0; + } + } + + /*--- Integration formulae from "Shape functions and points of + integration of the Résumé" by Josselin DELMAS (2013) ---*/ + + /*--- Triangle. Nodes of numerical integration at 1 point (order 1). ---*/ + + if (nNodes == 3) { + nGauss = 1; + Location[0][0] = 0.333333333333333; + Location[0][1] = 0.333333333333333; + Weight[0] = 0.5; + } + + /*--- Quadrilateral. Nodes of numerical integration at 4 points (order 2). ---*/ + + if (nNodes == 4) { + nGauss = 4; + Location[0][0] = -0.577350269189626; + Location[0][1] = -0.577350269189626; + Weight[0] = 1.0; + Location[1][0] = 0.577350269189626; + Location[1][1] = -0.577350269189626; + Weight[1] = 1.0; + Location[2][0] = 0.577350269189626; + Location[2][1] = 0.577350269189626; + Weight[2] = 1.0; + Location[3][0] = -0.577350269189626; + Location[3][1] = 0.577350269189626; + Weight[3] = 1.0; + } + + for (iGauss = 0; iGauss < nGauss; iGauss++) { + Xi = Location[iGauss][0]; + Eta = Location[iGauss][1]; + + if (nNodes == 3) Det = ShapeFunc_Triangle(Xi, Eta, CoordCorners, DShapeFunction); + if (nNodes == 4) Det = ShapeFunc_Quadrilateral(Xi, Eta, CoordCorners, DShapeFunction); + + /*--- Compute the B Matrix ---*/ + + for (iVar = 0; iVar < 3; iVar++) + for (jVar = 0; jVar < nNodes * nVar; jVar++) B_Matrix[iVar][jVar] = 0.0; + + for (iNode = 0; iNode < nNodes; iNode++) { + B_Matrix[0][0 + iNode * nVar] = DShapeFunction[iNode][0]; + B_Matrix[1][1 + iNode * nVar] = DShapeFunction[iNode][1]; + + B_Matrix[2][0 + iNode * nVar] = DShapeFunction[iNode][1]; + B_Matrix[2][1 + iNode * nVar] = DShapeFunction[iNode][0]; + } + + /*--- Impose a type of stiffness for each element ---*/ + + switch (config->GetDeform_Stiffness_Type()) { + case INVERSE_VOLUME: + E = 1.0 / ElemVolume; + break; + case SOLID_WALL_DISTANCE: + E = 1.0 / ElemDistance; + break; + case CONSTANT_STIFFNESS: + E = 1.0 / EPS; + break; + } + + Nu = config->GetDeform_Coeff(); + Mu = E / (2.0 * (1.0 + Nu)); + Lambda = Nu * E / ((1.0 + Nu) * (1.0 - 2.0 * Nu)); + + /*--- Compute the D Matrix (for plane strain and 3-D)---*/ + + D_Matrix[0][0] = Lambda + 2.0 * Mu; + D_Matrix[0][1] = Lambda; + D_Matrix[0][2] = 0.0; + D_Matrix[1][0] = Lambda; + D_Matrix[1][1] = Lambda + 2.0 * Mu; + D_Matrix[1][2] = 0.0; + D_Matrix[2][0] = 0.0; + D_Matrix[2][1] = 0.0; + D_Matrix[2][2] = Mu; + + /*--- Compute the BT.D Matrix ---*/ + + for (iVar = 0; iVar < nNodes * nVar; iVar++) { + for (jVar = 0; jVar < 3; jVar++) { + Aux_Matrix[iVar][jVar] = 0.0; + for (kVar = 0; kVar < 3; kVar++) Aux_Matrix[iVar][jVar] += B_Matrix[kVar][iVar] * D_Matrix[kVar][jVar]; + } + } + + /*--- Compute the BT.D.B Matrix (stiffness matrix), and add to the original + matrix using Gauss integration ---*/ + + for (iVar = 0; iVar < nNodes * nVar; iVar++) { + for (jVar = 0; jVar < nNodes * nVar; jVar++) { + for (kVar = 0; kVar < 3; kVar++) { + StiffMatrix_Elem[iVar][jVar] += Weight[iGauss] * Aux_Matrix[iVar][kVar] * B_Matrix[kVar][jVar] * fabs(Det); + } + } + } + } +} + + +void CLinearElasticity::SetFEA_StiffMatrix3D(CGeometry* geometry, CConfig* config, su2double** StiffMatrix_Elem, + unsigned long PointCorners[8], su2double CoordCorners[8][3], + unsigned short nNodes, su2double ElemVolume, su2double ElemDistance) { + su2double B_Matrix[6][24], + D_Matrix[6][6] = {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}, + Aux_Matrix[24][6]; + su2double Xi = 0.0, Eta = 0.0, Zeta = 0.0, Det = 0.0, Mu = 0.0, E = 0.0, Lambda = 0.0, Nu = 0.0; + unsigned short iNode, iVar, jVar, kVar, iGauss, nGauss = 0; + su2double DShapeFunction[8][4] = {{0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}}; + su2double Location[8][3], Weight[8]; + unsigned short nVar = geometry->GetnDim(); + + for (iVar = 0; iVar < nNodes * nVar; iVar++) { + for (jVar = 0; jVar < nNodes * nVar; jVar++) { + StiffMatrix_Elem[iVar][jVar] = 0.0; + } + } + + /*--- Integration formulae from "Shape functions and points of + integration of the Résumé" by Josselin Delmas (2013) ---*/ + + /*--- Tetrahedrons. Nodes of numerical integration at 1 point (order 1). ---*/ + + if (nNodes == 4) { + nGauss = 1; + Location[0][0] = 0.25; + Location[0][1] = 0.25; + Location[0][2] = 0.25; + Weight[0] = 0.166666666666666; + } + + /*--- Pyramids. Nodes numerical integration at 5 points. ---*/ + + if (nNodes == 5) { + nGauss = 5; + Location[0][0] = 0.5; + Location[0][1] = 0.0; + Location[0][2] = 0.1531754163448146; + Weight[0] = 0.133333333333333; + Location[1][0] = 0.0; + Location[1][1] = 0.5; + Location[1][2] = 0.1531754163448146; + Weight[1] = 0.133333333333333; + Location[2][0] = -0.5; + Location[2][1] = 0.0; + Location[2][2] = 0.1531754163448146; + Weight[2] = 0.133333333333333; + Location[3][0] = 0.0; + Location[3][1] = -0.5; + Location[3][2] = 0.1531754163448146; + Weight[3] = 0.133333333333333; + Location[4][0] = 0.0; + Location[4][1] = 0.0; + Location[4][2] = 0.6372983346207416; + Weight[4] = 0.133333333333333; + } + + /*--- Prism. Nodes of numerical integration at 6 points (order 3 in Xi, order 2 in Eta and Mu ). ---*/ + + if (nNodes == 6) { + nGauss = 6; + Location[0][0] = -0.577350269189626; + Location[0][1] = 0.166666666666667; + Location[0][2] = 0.166666666666667; + Weight[0] = 0.166666666666667; + Location[1][0] = -0.577350269189626; + Location[1][1] = 0.666666666666667; + Location[1][2] = 0.166666666666667; + Weight[1] = 0.166666666666667; + Location[2][0] = -0.577350269189626; + Location[2][1] = 0.166666666666667; + Location[2][2] = 0.666666666666667; + Weight[2] = 0.166666666666667; + Location[3][0] = 0.577350269189626; + Location[3][1] = 0.166666666666667; + Location[3][2] = 0.166666666666667; + Weight[3] = 0.166666666666667; + Location[4][0] = 0.577350269189626; + Location[4][1] = 0.666666666666667; + Location[4][2] = 0.166666666666667; + Weight[4] = 0.166666666666667; + Location[5][0] = 0.577350269189626; + Location[5][1] = 0.166666666666667; + Location[5][2] = 0.666666666666667; + Weight[5] = 0.166666666666667; + } + + /*--- Hexahedrons. Nodes of numerical integration at 6 points (order 3). ---*/ + + if (nNodes == 8) { + nGauss = 8; + Location[0][0] = -0.577350269189626; + Location[0][1] = -0.577350269189626; + Location[0][2] = -0.577350269189626; + Weight[0] = 1.0; + Location[1][0] = -0.577350269189626; + Location[1][1] = -0.577350269189626; + Location[1][2] = 0.577350269189626; + Weight[1] = 1.0; + Location[2][0] = -0.577350269189626; + Location[2][1] = 0.577350269189626; + Location[2][2] = -0.577350269189626; + Weight[2] = 1.0; + Location[3][0] = -0.577350269189626; + Location[3][1] = 0.577350269189626; + Location[3][2] = 0.577350269189626; + Weight[3] = 1.0; + Location[4][0] = 0.577350269189626; + Location[4][1] = -0.577350269189626; + Location[4][2] = -0.577350269189626; + Weight[4] = 1.0; + Location[5][0] = 0.577350269189626; + Location[5][1] = -0.577350269189626; + Location[5][2] = 0.577350269189626; + Weight[5] = 1.0; + Location[6][0] = 0.577350269189626; + Location[6][1] = 0.577350269189626; + Location[6][2] = -0.577350269189626; + Weight[6] = 1.0; + Location[7][0] = 0.577350269189626; + Location[7][1] = 0.577350269189626; + Location[7][2] = 0.577350269189626; + Weight[7] = 1.0; + } + + for (iGauss = 0; iGauss < nGauss; iGauss++) { + Xi = Location[iGauss][0]; + Eta = Location[iGauss][1]; + Zeta = Location[iGauss][2]; + + if (nNodes == 4) Det = ShapeFunc_Tetra(Xi, Eta, Zeta, CoordCorners, DShapeFunction); + if (nNodes == 5) Det = ShapeFunc_Pyram(Xi, Eta, Zeta, CoordCorners, DShapeFunction); + if (nNodes == 6) Det = ShapeFunc_Prism(Xi, Eta, Zeta, CoordCorners, DShapeFunction); + if (nNodes == 8) Det = ShapeFunc_Hexa(Xi, Eta, Zeta, CoordCorners, DShapeFunction); + + /*--- Compute the B Matrix ---*/ + + for (iVar = 0; iVar < 6; iVar++) + for (jVar = 0; jVar < nNodes * nVar; jVar++) B_Matrix[iVar][jVar] = 0.0; + + for (iNode = 0; iNode < nNodes; iNode++) { + B_Matrix[0][0 + iNode * nVar] = DShapeFunction[iNode][0]; + B_Matrix[1][1 + iNode * nVar] = DShapeFunction[iNode][1]; + B_Matrix[2][2 + iNode * nVar] = DShapeFunction[iNode][2]; + + B_Matrix[3][0 + iNode * nVar] = DShapeFunction[iNode][1]; + B_Matrix[3][1 + iNode * nVar] = DShapeFunction[iNode][0]; + + B_Matrix[4][1 + iNode * nVar] = DShapeFunction[iNode][2]; + B_Matrix[4][2 + iNode * nVar] = DShapeFunction[iNode][1]; + + B_Matrix[5][0 + iNode * nVar] = DShapeFunction[iNode][2]; + B_Matrix[5][2 + iNode * nVar] = DShapeFunction[iNode][0]; + } + + /*--- Impose a type of stiffness for each element ---*/ + + switch (config->GetDeform_Stiffness_Type()) { + case INVERSE_VOLUME: + E = 1.0 / ElemVolume; + break; + case SOLID_WALL_DISTANCE: + E = 1.0 / ElemDistance; + break; + case CONSTANT_STIFFNESS: + E = 1.0 / EPS; + break; + } + + Nu = config->GetDeform_Coeff(); + Mu = E / (2.0 * (1.0 + Nu)); + Lambda = Nu * E / ((1.0 + Nu) * (1.0 - 2.0 * Nu)); + + /*--- Compute the D Matrix (for plane strain and 3-D)---*/ + + D_Matrix[0][0] = Lambda + 2.0 * Mu; + D_Matrix[0][1] = Lambda; + D_Matrix[0][2] = Lambda; + D_Matrix[1][0] = Lambda; + D_Matrix[1][1] = Lambda + 2.0 * Mu; + D_Matrix[1][2] = Lambda; + D_Matrix[2][0] = Lambda; + D_Matrix[2][1] = Lambda; + D_Matrix[2][2] = Lambda + 2.0 * Mu; + D_Matrix[3][3] = Mu; + D_Matrix[4][4] = Mu; + D_Matrix[5][5] = Mu; + + /*--- Compute the BT.D Matrix ---*/ + + for (iVar = 0; iVar < nNodes * nVar; iVar++) { + for (jVar = 0; jVar < 6; jVar++) { + Aux_Matrix[iVar][jVar] = 0.0; + for (kVar = 0; kVar < 6; kVar++) Aux_Matrix[iVar][jVar] += B_Matrix[kVar][iVar] * D_Matrix[kVar][jVar]; + } + } + + /*--- Compute the BT.D.B Matrix (stiffness matrix), and add to the original + matrix using Gauss integration ---*/ + + for (iVar = 0; iVar < nNodes * nVar; iVar++) { + for (jVar = 0; jVar < nNodes * nVar; jVar++) { + for (kVar = 0; kVar < 6; kVar++) { + StiffMatrix_Elem[iVar][jVar] += Weight[iGauss] * Aux_Matrix[iVar][kVar] * B_Matrix[kVar][jVar] * fabs(Det); + } + } + } + } +} + + +void CLinearElasticity::AddFEA_StiffMatrix(CGeometry* geometry, su2double** StiffMatrix_Elem, + unsigned long PointCorners[8], unsigned short nNodes) { + unsigned short iVar, jVar, iDim, jDim; + + unsigned short nVar = geometry->GetnDim(); + + su2double** StiffMatrix_Node; + StiffMatrix_Node = new su2double*[nVar]; + for (iVar = 0; iVar < nVar; iVar++) StiffMatrix_Node[iVar] = new su2double[nVar]; + + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) StiffMatrix_Node[iVar][jVar] = 0.0; + + /*--- Transform the stiffness matrix for the hexahedral element into the + contributions for the individual nodes relative to each other. ---*/ + + for (iVar = 0; iVar < nNodes; iVar++) { + for (jVar = 0; jVar < nNodes; jVar++) { + for (iDim = 0; iDim < nVar; iDim++) { + for (jDim = 0; jDim < nVar; jDim++) { + StiffMatrix_Node[iDim][jDim] = StiffMatrix_Elem[(iVar * nVar) + iDim][(jVar * nVar) + jDim]; + } + } + + StiffMatrix.AddBlock(PointCorners[iVar], PointCorners[jVar], StiffMatrix_Node); + } + } + + /*--- Deallocate memory and exit ---*/ + + for (iVar = 0; iVar < nVar; iVar++) delete[] StiffMatrix_Node[iVar]; + delete[] StiffMatrix_Node; +} + +su2double CLinearElasticity::ShapeFunc_Triangle(su2double Xi, su2double Eta, su2double CoordCorners[8][3], + su2double DShapeFunction[8][4]) { + int i, j, k; + su2double c0, c1, xsj; + su2double xs[3][3], ad[3][3]; + + /*--- Shape functions ---*/ + + DShapeFunction[0][3] = Xi; + DShapeFunction[1][3] = Eta; + DShapeFunction[2][3] = 1 - Xi - Eta; + + /*--- dN/d xi, dN/d eta ---*/ + + DShapeFunction[0][0] = 1.0; + DShapeFunction[0][1] = 0.0; + DShapeFunction[1][0] = 0.0; + DShapeFunction[1][1] = 1.0; + DShapeFunction[2][0] = -1.0; + DShapeFunction[2][1] = -1.0; + + /*--- Jacobian transformation ---*/ + + for (i = 0; i < 2; i++) { + for (j = 0; j < 2; j++) { + xs[i][j] = 0.0; + for (k = 0; k < 3; k++) { + xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; + } + } + } + + /*--- Adjoint to Jacobian ---*/ + + ad[0][0] = xs[1][1]; + ad[0][1] = -xs[0][1]; + ad[1][0] = -xs[1][0]; + ad[1][1] = xs[0][0]; + + /*--- Determinant of Jacobian ---*/ + + xsj = ad[0][0] * ad[1][1] - ad[0][1] * ad[1][0]; + + /*--- Jacobian inverse ---*/ + + for (i = 0; i < 2; i++) { + for (j = 0; j < 2; j++) { + xs[i][j] = ad[i][j] / xsj; + } + } + + /*--- Derivatives with repect to global coordinates ---*/ + + for (k = 0; k < 3; k++) { + c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1]; // dN/dx + c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1]; // dN/dy + DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi + DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta + } + + return xsj; +} + +su2double CLinearElasticity::ShapeFunc_Quadrilateral(su2double Xi, su2double Eta, su2double CoordCorners[8][3], + su2double DShapeFunction[8][4]) { + int i, j, k; + su2double c0, c1, xsj; + su2double xs[3][3], ad[3][3]; + + /*--- Shape functions ---*/ + + DShapeFunction[0][3] = 0.25 * (1.0 - Xi) * (1.0 - Eta); + DShapeFunction[1][3] = 0.25 * (1.0 + Xi) * (1.0 - Eta); + DShapeFunction[2][3] = 0.25 * (1.0 + Xi) * (1.0 + Eta); + DShapeFunction[3][3] = 0.25 * (1.0 - Xi) * (1.0 + Eta); + + /*--- dN/d xi, dN/d eta ---*/ + + DShapeFunction[0][0] = -0.25 * (1.0 - Eta); + DShapeFunction[0][1] = -0.25 * (1.0 - Xi); + DShapeFunction[1][0] = 0.25 * (1.0 - Eta); + DShapeFunction[1][1] = -0.25 * (1.0 + Xi); + DShapeFunction[2][0] = 0.25 * (1.0 + Eta); + DShapeFunction[2][1] = 0.25 * (1.0 + Xi); + DShapeFunction[3][0] = -0.25 * (1.0 + Eta); + DShapeFunction[3][1] = 0.25 * (1.0 - Xi); + + /*--- Jacobian transformation ---*/ + + for (i = 0; i < 2; i++) { + for (j = 0; j < 2; j++) { + xs[i][j] = 0.0; + for (k = 0; k < 4; k++) { + xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; + } + } + } + + /*--- Adjoint to Jacobian ---*/ + + ad[0][0] = xs[1][1]; + ad[0][1] = -xs[0][1]; + ad[1][0] = -xs[1][0]; + ad[1][1] = xs[0][0]; + + /*--- Determinant of Jacobian ---*/ + + xsj = ad[0][0] * ad[1][1] - ad[0][1] * ad[1][0]; + + /*--- Jacobian inverse ---*/ + + for (i = 0; i < 2; i++) { + for (j = 0; j < 2; j++) { + xs[i][j] = ad[i][j] / xsj; + } + } + + /*--- Derivatives with repect to global coordinates ---*/ + + for (k = 0; k < 4; k++) { + c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1]; // dN/dx + c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1]; // dN/dy + DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi + DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta + } + + return xsj; +} + +su2double CLinearElasticity::ShapeFunc_Tetra(su2double Xi, su2double Eta, su2double Zeta, + su2double CoordCorners[8][3], su2double DShapeFunction[8][4]) { + int i, j, k; + su2double c0, c1, c2, xsj; + su2double xs[3][3], ad[3][3]; + + /*--- Shape functions ---*/ + + DShapeFunction[0][3] = Xi; + DShapeFunction[1][3] = Zeta; + DShapeFunction[2][3] = 1.0 - Xi - Eta - Zeta; + DShapeFunction[3][3] = Eta; + + /*--- dN/d xi, dN/d eta, dN/d zeta ---*/ + + DShapeFunction[0][0] = 1.0; + DShapeFunction[0][1] = 0.0; + DShapeFunction[0][2] = 0.0; + DShapeFunction[1][0] = 0.0; + DShapeFunction[1][1] = 0.0; + DShapeFunction[1][2] = 1.0; + DShapeFunction[2][0] = -1.0; + DShapeFunction[2][1] = -1.0; + DShapeFunction[2][2] = -1.0; + DShapeFunction[3][0] = 0.0; + DShapeFunction[3][1] = 1.0; + DShapeFunction[3][2] = 0.0; + + /*--- Jacobian transformation ---*/ + + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + xs[i][j] = 0.0; + for (k = 0; k < 4; k++) { + xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; + } + } + } + + /*--- Adjoint to Jacobian ---*/ + + ad[0][0] = xs[1][1] * xs[2][2] - xs[1][2] * xs[2][1]; + ad[0][1] = xs[0][2] * xs[2][1] - xs[0][1] * xs[2][2]; + ad[0][2] = xs[0][1] * xs[1][2] - xs[0][2] * xs[1][1]; + ad[1][0] = xs[1][2] * xs[2][0] - xs[1][0] * xs[2][2]; + ad[1][1] = xs[0][0] * xs[2][2] - xs[0][2] * xs[2][0]; + ad[1][2] = xs[0][2] * xs[1][0] - xs[0][0] * xs[1][2]; + ad[2][0] = xs[1][0] * xs[2][1] - xs[1][1] * xs[2][0]; + ad[2][1] = xs[0][1] * xs[2][0] - xs[0][0] * xs[2][1]; + ad[2][2] = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; + + /*--- Determinant of Jacobian ---*/ + + xsj = xs[0][0] * ad[0][0] + xs[0][1] * ad[1][0] + xs[0][2] * ad[2][0]; + + /*--- Jacobian inverse ---*/ + + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + xs[i][j] = ad[i][j] / xsj; + } + } + + /*--- Derivatives with repect to global coordinates ---*/ + + for (k = 0; k < 4; k++) { + c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1] + xs[0][2] * DShapeFunction[k][2]; // dN/dx + c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1] + xs[1][2] * DShapeFunction[k][2]; // dN/dy + c2 = xs[2][0] * DShapeFunction[k][0] + xs[2][1] * DShapeFunction[k][1] + xs[2][2] * DShapeFunction[k][2]; // dN/dz + DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi + DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta + DShapeFunction[k][2] = c2; // store dN/dz instead of dN/d zeta + } + + return xsj; +} + +su2double CLinearElasticity::ShapeFunc_Pyram(su2double Xi, su2double Eta, su2double Zeta, + su2double CoordCorners[8][3], su2double DShapeFunction[8][4]) { + int i, j, k; + su2double c0, c1, c2, xsj; + su2double xs[3][3], ad[3][3]; + + /*--- Shape functions ---*/ + + DShapeFunction[0][3] = 0.25 * (-Xi + Eta + Zeta - 1.0) * (-Xi - Eta + Zeta - 1.0) / (1.0 - Zeta); + DShapeFunction[1][3] = 0.25 * (-Xi - Eta + Zeta - 1.0) * (Xi - Eta + Zeta - 1.0) / (1.0 - Zeta); + DShapeFunction[2][3] = 0.25 * (Xi + Eta + Zeta - 1.0) * (Xi - Eta + Zeta - 1.0) / (1.0 - Zeta); + DShapeFunction[3][3] = 0.25 * (Xi + Eta + Zeta - 1.0) * (-Xi + Eta + Zeta - 1.0) / (1.0 - Zeta); + DShapeFunction[4][3] = Zeta; + + /*--- dN/d xi ---*/ + + DShapeFunction[0][0] = 0.5 * (Zeta - Xi - 1.0) / (Zeta - 1.0); + DShapeFunction[1][0] = 0.5 * Xi / (Zeta - 1.0); + DShapeFunction[2][0] = 0.5 * (1.0 - Zeta - Xi) / (Zeta - 1.0); + DShapeFunction[3][0] = DShapeFunction[1][0]; + DShapeFunction[4][0] = 0.0; + + /*--- dN/d eta ---*/ + + DShapeFunction[0][1] = 0.5 * Eta / (Zeta - 1.0); + DShapeFunction[1][1] = 0.5 * (Zeta - Eta - 1.0) / (Zeta - 1.0); + DShapeFunction[2][1] = DShapeFunction[0][1]; + DShapeFunction[3][1] = 0.5 * (1.0 - Zeta - Eta) / (Zeta - 1.0); + DShapeFunction[4][1] = 0.0; + + /*--- dN/d zeta ---*/ + + DShapeFunction[0][2] = 0.25 * (-1.0 + 2.0 * Zeta - Zeta * Zeta - Eta * Eta + Xi * Xi) / ((1.0 - Zeta) * (1.0 - Zeta)); + DShapeFunction[1][2] = 0.25 * (-1.0 + 2.0 * Zeta - Zeta * Zeta + Eta * Eta - Xi * Xi) / ((1.0 - Zeta) * (1.0 - Zeta)); + DShapeFunction[2][2] = DShapeFunction[0][2]; + DShapeFunction[3][2] = DShapeFunction[1][2]; + DShapeFunction[4][2] = 1.0; + + /*--- Jacobian transformation ---*/ + + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + xs[i][j] = 0.0; + for (k = 0; k < 5; k++) { + xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; + } + } + } + + /*--- Adjoint to Jacobian ---*/ + + ad[0][0] = xs[1][1] * xs[2][2] - xs[1][2] * xs[2][1]; + ad[0][1] = xs[0][2] * xs[2][1] - xs[0][1] * xs[2][2]; + ad[0][2] = xs[0][1] * xs[1][2] - xs[0][2] * xs[1][1]; + ad[1][0] = xs[1][2] * xs[2][0] - xs[1][0] * xs[2][2]; + ad[1][1] = xs[0][0] * xs[2][2] - xs[0][2] * xs[2][0]; + ad[1][2] = xs[0][2] * xs[1][0] - xs[0][0] * xs[1][2]; + ad[2][0] = xs[1][0] * xs[2][1] - xs[1][1] * xs[2][0]; + ad[2][1] = xs[0][1] * xs[2][0] - xs[0][0] * xs[2][1]; + ad[2][2] = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; + + /*--- Determinant of Jacobian ---*/ + + xsj = xs[0][0] * ad[0][0] + xs[0][1] * ad[1][0] + xs[0][2] * ad[2][0]; + + /*--- Jacobian inverse ---*/ + + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + xs[i][j] = ad[i][j] / xsj; + } + } + + /*--- Derivatives with repect to global coordinates ---*/ + + for (k = 0; k < 5; k++) { + c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1] + xs[0][2] * DShapeFunction[k][2]; // dN/dx + c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1] + xs[1][2] * DShapeFunction[k][2]; // dN/dy + c2 = xs[2][0] * DShapeFunction[k][0] + xs[2][1] * DShapeFunction[k][1] + xs[2][2] * DShapeFunction[k][2]; // dN/dz + DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi + DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta + DShapeFunction[k][2] = c2; // store dN/dz instead of dN/d zeta + } + + return xsj; +} + +su2double CLinearElasticity::ShapeFunc_Prism(su2double Xi, su2double Eta, su2double Zeta, + su2double CoordCorners[8][3], su2double DShapeFunction[8][4]) { + int i, j, k; + su2double c0, c1, c2, xsj; + su2double xs[3][3], ad[3][3]; + + /*--- Shape functions ---*/ + + DShapeFunction[0][3] = 0.5 * Eta * (1.0 - Xi); + DShapeFunction[1][3] = 0.5 * Zeta * (1.0 - Xi); + DShapeFunction[2][3] = 0.5 * (1.0 - Eta - Zeta) * (1.0 - Xi); + DShapeFunction[3][3] = 0.5 * Eta * (Xi + 1.0); + DShapeFunction[4][3] = 0.5 * Zeta * (Xi + 1.0); + DShapeFunction[5][3] = 0.5 * (1.0 - Eta - Zeta) * (Xi + 1.0); + + /*--- dN/d Xi, dN/d Eta, dN/d Zeta ---*/ + + DShapeFunction[0][0] = -0.5 * Eta; + DShapeFunction[0][1] = 0.5 * (1.0 - Xi); + DShapeFunction[0][2] = 0.0; + DShapeFunction[1][0] = -0.5 * Zeta; + DShapeFunction[1][1] = 0.0; + DShapeFunction[1][2] = 0.5 * (1.0 - Xi); + DShapeFunction[2][0] = -0.5 * (1.0 - Eta - Zeta); + DShapeFunction[2][1] = -0.5 * (1.0 - Xi); + DShapeFunction[2][2] = -0.5 * (1.0 - Xi); + DShapeFunction[3][0] = 0.5 * Eta; + DShapeFunction[3][1] = 0.5 * (Xi + 1.0); + DShapeFunction[3][2] = 0.0; + DShapeFunction[4][0] = 0.5 * Zeta; + DShapeFunction[4][1] = 0.0; + DShapeFunction[4][2] = 0.5 * (Xi + 1.0); + DShapeFunction[5][0] = 0.5 * (1.0 - Eta - Zeta); + DShapeFunction[5][1] = -0.5 * (Xi + 1.0); + DShapeFunction[5][2] = -0.5 * (Xi + 1.0); + + /*--- Jacobian transformation ---*/ + + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + xs[i][j] = 0.0; + for (k = 0; k < 6; k++) { + xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; + } + } + } + + /*--- Adjoint to Jacobian ---*/ + + ad[0][0] = xs[1][1] * xs[2][2] - xs[1][2] * xs[2][1]; + ad[0][1] = xs[0][2] * xs[2][1] - xs[0][1] * xs[2][2]; + ad[0][2] = xs[0][1] * xs[1][2] - xs[0][2] * xs[1][1]; + ad[1][0] = xs[1][2] * xs[2][0] - xs[1][0] * xs[2][2]; + ad[1][1] = xs[0][0] * xs[2][2] - xs[0][2] * xs[2][0]; + ad[1][2] = xs[0][2] * xs[1][0] - xs[0][0] * xs[1][2]; + ad[2][0] = xs[1][0] * xs[2][1] - xs[1][1] * xs[2][0]; + ad[2][1] = xs[0][1] * xs[2][0] - xs[0][0] * xs[2][1]; + ad[2][2] = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; + + /*--- Determinant of Jacobian ---*/ + + xsj = xs[0][0] * ad[0][0] + xs[0][1] * ad[1][0] + xs[0][2] * ad[2][0]; + + /*--- Jacobian inverse ---*/ + + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + xs[i][j] = ad[i][j] / xsj; + } + } + + /*--- Derivatives with repect to global coordinates ---*/ + + for (k = 0; k < 6; k++) { + c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1] + xs[0][2] * DShapeFunction[k][2]; // dN/dx + c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1] + xs[1][2] * DShapeFunction[k][2]; // dN/dy + c2 = xs[2][0] * DShapeFunction[k][0] + xs[2][1] * DShapeFunction[k][1] + xs[2][2] * DShapeFunction[k][2]; // dN/dz + DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi + DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta + DShapeFunction[k][2] = c2; // store dN/dz instead of dN/d zeta + } + + return xsj; +} + +su2double CLinearElasticity::ShapeFunc_Hexa(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3], + su2double DShapeFunction[8][4]) { + int i, j, k; + su2double c0, c1, c2, xsj; + su2double xs[3][3], ad[3][3]; + + /*--- Shape functions ---*/ + + DShapeFunction[0][3] = 0.125 * (1.0 - Xi) * (1.0 - Eta) * (1.0 - Zeta); + DShapeFunction[1][3] = 0.125 * (1.0 + Xi) * (1.0 - Eta) * (1.0 - Zeta); + DShapeFunction[2][3] = 0.125 * (1.0 + Xi) * (1.0 + Eta) * (1.0 - Zeta); + DShapeFunction[3][3] = 0.125 * (1.0 - Xi) * (1.0 + Eta) * (1.0 - Zeta); + DShapeFunction[4][3] = 0.125 * (1.0 - Xi) * (1.0 - Eta) * (1.0 + Zeta); + DShapeFunction[5][3] = 0.125 * (1.0 + Xi) * (1.0 - Eta) * (1.0 + Zeta); + DShapeFunction[6][3] = 0.125 * (1.0 + Xi) * (1.0 + Eta) * (1.0 + Zeta); + DShapeFunction[7][3] = 0.125 * (1.0 - Xi) * (1.0 + Eta) * (1.0 + Zeta); + + /*--- dN/d xi ---*/ + + DShapeFunction[0][0] = -0.125 * (1.0 - Eta) * (1.0 - Zeta); + DShapeFunction[1][0] = 0.125 * (1.0 - Eta) * (1.0 - Zeta); + DShapeFunction[2][0] = 0.125 * (1.0 + Eta) * (1.0 - Zeta); + DShapeFunction[3][0] = -0.125 * (1.0 + Eta) * (1.0 - Zeta); + DShapeFunction[4][0] = -0.125 * (1.0 - Eta) * (1.0 + Zeta); + DShapeFunction[5][0] = 0.125 * (1.0 - Eta) * (1.0 + Zeta); + DShapeFunction[6][0] = 0.125 * (1.0 + Eta) * (1.0 + Zeta); + DShapeFunction[7][0] = -0.125 * (1.0 + Eta) * (1.0 + Zeta); + + /*--- dN/d eta ---*/ + + DShapeFunction[0][1] = -0.125 * (1.0 - Xi) * (1.0 - Zeta); + DShapeFunction[1][1] = -0.125 * (1.0 + Xi) * (1.0 - Zeta); + DShapeFunction[2][1] = 0.125 * (1.0 + Xi) * (1.0 - Zeta); + DShapeFunction[3][1] = 0.125 * (1.0 - Xi) * (1.0 - Zeta); + DShapeFunction[4][1] = -0.125 * (1.0 - Xi) * (1.0 + Zeta); + DShapeFunction[5][1] = -0.125 * (1.0 + Xi) * (1.0 + Zeta); + DShapeFunction[6][1] = 0.125 * (1.0 + Xi) * (1.0 + Zeta); + DShapeFunction[7][1] = 0.125 * (1.0 - Xi) * (1.0 + Zeta); + + /*--- dN/d zeta ---*/ + + DShapeFunction[0][2] = -0.125 * (1.0 - Xi) * (1.0 - Eta); + DShapeFunction[1][2] = -0.125 * (1.0 + Xi) * (1.0 - Eta); + DShapeFunction[2][2] = -0.125 * (1.0 + Xi) * (1.0 + Eta); + DShapeFunction[3][2] = -0.125 * (1.0 - Xi) * (1.0 + Eta); + DShapeFunction[4][2] = 0.125 * (1.0 - Xi) * (1.0 - Eta); + DShapeFunction[5][2] = 0.125 * (1.0 + Xi) * (1.0 - Eta); + DShapeFunction[6][2] = 0.125 * (1.0 + Xi) * (1.0 + Eta); + DShapeFunction[7][2] = 0.125 * (1.0 - Xi) * (1.0 + Eta); + + /*--- Jacobian transformation ---*/ + + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + xs[i][j] = 0.0; + for (k = 0; k < 8; k++) { + xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; + } + } + } + + /*--- Adjoint to Jacobian ---*/ + + ad[0][0] = xs[1][1] * xs[2][2] - xs[1][2] * xs[2][1]; + ad[0][1] = xs[0][2] * xs[2][1] - xs[0][1] * xs[2][2]; + ad[0][2] = xs[0][1] * xs[1][2] - xs[0][2] * xs[1][1]; + ad[1][0] = xs[1][2] * xs[2][0] - xs[1][0] * xs[2][2]; + ad[1][1] = xs[0][0] * xs[2][2] - xs[0][2] * xs[2][0]; + ad[1][2] = xs[0][2] * xs[1][0] - xs[0][0] * xs[1][2]; + ad[2][0] = xs[1][0] * xs[2][1] - xs[1][1] * xs[2][0]; + ad[2][1] = xs[0][1] * xs[2][0] - xs[0][0] * xs[2][1]; + ad[2][2] = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; + + /*--- Determinant of Jacobian ---*/ + + xsj = xs[0][0] * ad[0][0] + xs[0][1] * ad[1][0] + xs[0][2] * ad[2][0]; + + /*--- Jacobian inverse ---*/ + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + xs[i][j] = ad[i][j] / xsj; + } + } + + /*--- Derivatives with repect to global coordinates ---*/ + + for (k = 0; k < 8; k++) { + c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1] + xs[0][2] * DShapeFunction[k][2]; // dN/dx + c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1] + xs[1][2] * DShapeFunction[k][2]; // dN/dy + c2 = xs[2][0] * DShapeFunction[k][0] + xs[2][1] * DShapeFunction[k][1] + xs[2][2] * DShapeFunction[k][2]; // dN/dz + DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi + DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta + DShapeFunction[k][2] = c2; // store dN/dz instead of dN/d zeta + } + + return xsj; +} + +void CLinearElasticity::SetDomainDisplacements(CGeometry* geometry, CConfig* config) { + unsigned short iDim, nDim = geometry->GetnDim(); + unsigned long iPoint, total_index; + + if (config->GetHold_GridFixed()) { + auto MinCoordValues = config->GetHold_GridFixed_Coord(); + auto MaxCoordValues = &config->GetHold_GridFixed_Coord()[3]; + + /*--- Set to zero displacements of all the points that are not going to be moved + except the surfaces ---*/ + + for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + auto Coord = geometry->nodes->GetCoord(iPoint); + for (iDim = 0; iDim < nDim; iDim++) { + if ((Coord[iDim] < MinCoordValues[iDim]) || (Coord[iDim] > MaxCoordValues[iDim])) { + total_index = iPoint * nDim + iDim; + LinSysRes[total_index] = 0.0; + LinSysSol[total_index] = 0.0; + StiffMatrix.DeleteValsRowi(total_index); + } + } + } + } + + /*--- Don't move the volume grid outside the limits based + on the distance to the solid surface ---*/ + + if (config->GetDeform_Limit() < 1E6) { + for (iPoint = 0; iPoint < nPoint; iPoint++) { + if (geometry->nodes->GetWall_Distance(iPoint) >= config->GetDeform_Limit()) { + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + LinSysRes[total_index] = 0.0; + LinSysSol[total_index] = 0.0; + StiffMatrix.DeleteValsRowi(total_index); + } + } + } + } +} + + +void CLinearElasticity::SetBoundaryDisplacements(CGeometry* geometry, CConfig* config) { + unsigned short iDim, nDim = geometry->GetnDim(), iMarker, axis = 0; + unsigned long iPoint, total_index, iVertex; + su2double *VarCoord, MeanCoord[3] = {0.0, 0.0, 0.0}, VarIncrement = 1.0; + + /*--- Get the SU2 module. SU2_CFD will use this routine for dynamically + deforming meshes (MARKER_MOVING), while SU2_DEF will use it for deforming + meshes after imposing design variable surface deformations (DV_MARKER). ---*/ + + SU2_COMPONENT Kind_SU2 = config->GetKind_SU2(); + + /*--- If requested (no by default) impose the surface deflections in + increments and solve the grid deformation equations iteratively with + successive small deformations. ---*/ + + VarIncrement = 1.0 / ((su2double)config->GetGridDef_Nonlinear_Iter()); + + /*--- As initialization, set to zero displacements of all the surfaces except the symmetry + plane (which is treated specially, see below), internal and the send-receive boundaries ---*/ + + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if (((config->GetMarker_All_KindBC(iMarker) != SYMMETRY_PLANE) && + (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) && + (config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY))) { + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + LinSysRes[total_index] = 0.0; + LinSysSol[total_index] = 0.0; + StiffMatrix.DeleteValsRowi(total_index); + } + } + } + } + + /*--- Set the known displacements, note that some points of the moving surfaces + could be on on the symmetry plane, we should specify DeleteValsRowi again (just in case) ---*/ + + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if (((config->GetMarker_All_Moving(iMarker) == YES) && (Kind_SU2 == SU2_COMPONENT::SU2_CFD)) || + ((config->GetMarker_All_DV(iMarker) == YES) && (Kind_SU2 == SU2_COMPONENT::SU2_DEF)) || + ((config->GetDirectDiff() == D_DESIGN) && (Kind_SU2 == SU2_COMPONENT::SU2_CFD) && + (config->GetMarker_All_DV(iMarker) == YES)) || + ((config->GetMarker_All_DV(iMarker) == YES) && (Kind_SU2 == SU2_COMPONENT::SU2_DOT))) { + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + + VarCoord = geometry->vertex[iMarker][iVertex]->GetVarCoord(); + + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + LinSysRes[total_index] = SU2_TYPE::GetValue(VarCoord[iDim] * VarIncrement); + LinSysSol[total_index] = SU2_TYPE::GetValue(VarCoord[iDim] * VarIncrement); + StiffMatrix.DeleteValsRowi(total_index); + } + } + } + } + + /*--- Set to zero displacements of the normal component for the symmetry plane condition ---*/ + + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if ((config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE)) { + su2double* Coord_0 = nullptr; + + for (iDim = 0; iDim < nDim; iDim++) MeanCoord[iDim] = 0.0; + + /*--- Store the coord of the first point to help identify the axis. ---*/ + + iPoint = geometry->vertex[iMarker][0]->GetNode(); + Coord_0 = geometry->nodes->GetCoord(iPoint); + + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + VarCoord = geometry->nodes->GetCoord(iPoint); + for (iDim = 0; iDim < nDim; iDim++) + MeanCoord[iDim] += (VarCoord[iDim] - Coord_0[iDim]) * (VarCoord[iDim] - Coord_0[iDim]); + } + for (iDim = 0; iDim < nDim; iDim++) MeanCoord[iDim] = sqrt(MeanCoord[iDim]); + if (nDim == 3) { + if ((MeanCoord[0] <= MeanCoord[1]) && (MeanCoord[0] <= MeanCoord[2])) axis = 0; + if ((MeanCoord[1] <= MeanCoord[0]) && (MeanCoord[1] <= MeanCoord[2])) axis = 1; + if ((MeanCoord[2] <= MeanCoord[0]) && (MeanCoord[2] <= MeanCoord[1])) axis = 2; + } else { + if ((MeanCoord[0] <= MeanCoord[1])) axis = 0; + if ((MeanCoord[1] <= MeanCoord[0])) axis = 1; + } + + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + total_index = iPoint * nDim + axis; + LinSysRes[total_index] = 0.0; + LinSysSol[total_index] = 0.0; + StiffMatrix.DeleteValsRowi(total_index); + } + } + } + + /*--- Don't move the nearfield plane ---*/ + + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == NEARFIELD_BOUNDARY) { + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + LinSysRes[total_index] = 0.0; + LinSysSol[total_index] = 0.0; + StiffMatrix.DeleteValsRowi(total_index); + } + } + } + } + + /*--- Move the FSI interfaces ---*/ + + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if ((config->GetMarker_All_ZoneInterface(iMarker) == YES) && (Kind_SU2 == SU2_COMPONENT::SU2_CFD)) { + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + VarCoord = geometry->vertex[iMarker][iVertex]->GetVarCoord(); + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + LinSysRes[total_index] = SU2_TYPE::GetValue(VarCoord[iDim] * VarIncrement); + LinSysSol[total_index] = SU2_TYPE::GetValue(VarCoord[iDim] * VarIncrement); + StiffMatrix.DeleteValsRowi(total_index); + } + } + } + } +} + +void CLinearElasticity::SetBoundaryDerivatives(CGeometry* geometry, CConfig* config, + bool ForwardProjectionDerivative) { + unsigned short iDim, iMarker; + unsigned long iPoint, total_index, iVertex; + + su2double* VarCoord; + SU2_COMPONENT Kind_SU2 = config->GetKind_SU2(); + if ((config->GetDirectDiff() == D_DESIGN) && (Kind_SU2 == SU2_COMPONENT::SU2_CFD)) { + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if ((config->GetMarker_All_DV(iMarker) == YES)) { + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + VarCoord = geometry->vertex[iMarker][iVertex]->GetVarCoord(); + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + LinSysRes[total_index] = SU2_TYPE::GetDerivative(VarCoord[iDim]); + LinSysSol[total_index] = SU2_TYPE::GetDerivative(VarCoord[iDim]); + } + } + } + } + if (LinSysRes.norm() == 0.0) cout << "Warning: Derivatives are zero!" << endl; + } else if ((Kind_SU2 == SU2_COMPONENT::SU2_DOT) && !ForwardProjectionDerivative) { + for (iPoint = 0; iPoint < nPoint; iPoint++) { + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + LinSysRes[total_index] = SU2_TYPE::GetValue(geometry->GetSensitivity(iPoint, iDim)); + LinSysSol[total_index] = SU2_TYPE::GetValue(geometry->GetSensitivity(iPoint, iDim)); + } + } + } else if (config->GetSmoothGradient() && ForwardProjectionDerivative) { + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if ((config->GetMarker_All_DV(iMarker) == YES)) { + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + for (iDim = 0; iDim < nDim; iDim++) { + total_index = iPoint * nDim + iDim; + LinSysRes[total_index] = SU2_TYPE::GetValue(geometry->GetSensitivity(iPoint, iDim)); + LinSysSol[total_index] = SU2_TYPE::GetValue(geometry->GetSensitivity(iPoint, iDim)); + } + } + } + } + if (LinSysRes.norm() == 0.0) cout << "Warning: Derivatives are zero!" << endl; + } +} diff --git a/Common/src/grid_movement/CRadialBasisFunctionInterpolation.cpp b/Common/src/grid_movement/CRadialBasisFunctionInterpolation.cpp new file mode 100644 index 000000000000..6ec015b6e273 --- /dev/null +++ b/Common/src/grid_movement/CRadialBasisFunctionInterpolation.cpp @@ -0,0 +1,715 @@ +/*! + * \file CRadialBasisFunctionInterpolation.cpp + * \brief Subroutines for moving mesh volume elements using Radial Basis Function interpolation. + * \author F. van Steen + * \version 8.0.1 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/grid_movement/CRadialBasisFunctionInterpolation.hpp" +#include "../../include/interface_interpolation/CRadialBasisFunction.hpp" +#include "../../include/toolboxes/geometry_toolbox.hpp" +#include "../../include/adt/CADTPointsOnlyClass.hpp" + + +CRadialBasisFunctionInterpolation::CRadialBasisFunctionInterpolation(CGeometry* geometry, CConfig* config) : CVolumetricMovement(geometry) {} + +CRadialBasisFunctionInterpolation::~CRadialBasisFunctionInterpolation() = default; + +void CRadialBasisFunctionInterpolation::SetVolume_Deformation(CGeometry* geometry, CConfig* config, bool UpdateGeo, bool Derivative, + bool ForwardProjectionDerivative){ + + + /*--- Retrieve type of RBF and its support radius ---*/ + + const auto kindRBF = config->GetKindRadialBasisFunction(); + const su2double radius = config->GetRadialBasisFunctionParameter(); + + su2double MinVolume, MaxVolume; + + /*--- Retrieving number of deformation steps and screen output from config ---*/ + + const auto Nonlinear_Iter = config->GetGridDef_Nonlinear_Iter(); + auto Screen_Output = config->GetDeform_Output(); + + /*--- Disable the screen output if we're running SU2_CFD ---*/ + + if (config->GetKind_SU2() == SU2_COMPONENT::SU2_CFD && !Derivative) Screen_Output = false; + if (config->GetSmoothGradient()) Screen_Output = true; + + + /*--- Determining the boundary and internal nodes. Setting the control nodes. ---*/ + SetBoundNodes(geometry, config); + + vector internalNodes; + SetInternalNodes(geometry, config, internalNodes); + + SetCtrlNodes(config); + + /*--- Looping over the number of deformation iterations ---*/ + for (auto iNonlinear_Iter = 0ul; iNonlinear_Iter < Nonlinear_Iter; iNonlinear_Iter++) { + + /*--- Compute min volume in the entire mesh. ---*/ + + ComputeDeforming_Element_Volume(geometry, MinVolume, MaxVolume, Screen_Output); + if (rank == MASTER_NODE && Screen_Output) + cout << "Min. volume: " << MinVolume << ", max. volume: " << MaxVolume << "." << endl; + + + /*--- Solving the RBF system, resulting in the interpolation coefficients ---*/ + SolveRBF_System(geometry, config, kindRBF, radius); + + /*--- Updating the coordinates of the grid ---*/ + UpdateGridCoord(geometry, config, kindRBF, radius, internalNodes); + + if(UpdateGeo){ + UpdateDualGrid(geometry, config); + } + + /*--- Check for failed deformation (negative volumes). ---*/ + + ComputeDeforming_Element_Volume(geometry, MinVolume, MaxVolume, Screen_Output); + + /*--- Calculate amount of nonconvex elements ---*/ + + ComputenNonconvexElements(geometry, Screen_Output); + + if (rank == MASTER_NODE && Screen_Output) { + cout << "Non-linear iter.: " << iNonlinear_Iter + 1 << "/" << Nonlinear_Iter << ". "; + if (nDim == 2) + cout << "Min. area: " << MinVolume << "." << endl; + else + cout << "Min. volume: " << MinVolume << "." << endl; + } + } +} + +void CRadialBasisFunctionInterpolation::SolveRBF_System(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius){ + + /*--- In case of data reduction an iterative greedy algorithm is applied + to perform the interpolation with a reduced set of control nodes. + Otherwise with a full set of control nodes. ---*/ + + if(config->GetRBF_DataReduction()){ + + /*--- Error tolerance for the data reduction tolerance ---*/ + const su2double dataReductionTolerance = config->GetRBF_DataRedTolerance(); + + /*--- Local maximum error node and corresponding maximum error ---*/ + unsigned long maxErrorNodeLocal; + su2double maxErrorLocal{0}; + + /*--- Obtaining the initial maximum error nodes, which are found based on the maximum applied deformation. */ + if(ControlNodes->empty()){ + GetInitMaxErrorNode(geometry, config, maxErrorNodeLocal, maxErrorLocal); + SU2_MPI::Allreduce(&maxErrorLocal, &MaxErrorGlobal, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); + } + + /*--- Number of greedy iterations. ---*/ + unsigned short greedyIter = 0; + + /*--- While the maximum error is above the tolerance, data reduction algorithm is continued. ---*/ + while(MaxErrorGlobal > dataReductionTolerance || greedyIter == 0){ + + /*--- In case of a nonzero local error, control nodes are added ---*/ + if(maxErrorLocal> 0){ + AddControlNode(maxErrorNodeLocal); + } + + /*--- Obtaining the global number of control nodes. ---*/ + Get_nCtrlNodesGlobal(); + + /*--- Obtaining the interpolation coefficients. ---*/ + GetInterpCoeffs(geometry, config, type, radius); + + /*--- Determining the interpolation error, of the non-control boundary nodes. ---*/ + GetInterpError(geometry, config, type, radius, maxErrorNodeLocal, maxErrorLocal); + SU2_MPI::Allreduce(&maxErrorLocal, &MaxErrorGlobal, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); + + if(rank == MASTER_NODE) cout << "Greedy iteration: " << greedyIter << ". Max error: " << MaxErrorGlobal << ". Global nr. of ctrl nodes: " << nCtrlNodesGlobal << "\n" << endl; + greedyIter++; + + } + }else{ + + /*--- Obtaining the interpolation coefficients. ---*/ + GetInterpCoeffs(geometry, config, type, radius); + } +} + +void CRadialBasisFunctionInterpolation::GetInterpCoeffs(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius){ + + /*--- Obtaining the control nodes coordinates and distributing over all processes. ---*/ + SetCtrlNodeCoords(geometry); + + /*--- Obtaining the deformation of the control nodes. ---*/ + SetDeformation(geometry, config); + + /*--- Computation of the (inverse) interpolation matrix. ---*/ + su2passivematrix invInterpMat; + ComputeInterpolationMatrix(geometry, type, radius, invInterpMat); + + /*--- Obtaining the interpolation coefficients. ---*/ + ComputeInterpCoeffs(invInterpMat); +} + + +void CRadialBasisFunctionInterpolation::SetBoundNodes(CGeometry* geometry, CConfig* config){ + + /*--- Storing of the local node, marker and vertex information of the boundary nodes ---*/ + + /*--- Looping over the markers ---*/ + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { + + /*--- Checking if not internal or send/receive marker ---*/ + if (!config->GetMarker_All_Deform_Mesh_Internal(iMarker) && !config->GetMarker_All_SendRecv(iMarker)) { + + /*--- Looping over the vertices of marker ---*/ + for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) { + + /*--- Node in consideration ---*/ + auto iNode = geometry->vertex[iMarker][iVertex]->GetNode(); + + /*--- Check whether node is part of the subdomain and not shared with a receiving marker (for parallel computation)*/ + if (geometry->nodes->GetDomain(iNode)) { + BoundNodes.push_back(new CRadialBasisFunctionNode(iNode, iMarker, iVertex)); + } + } + } + } + + /*--- Sorting of the boundary nodes based on their index ---*/ + sort(BoundNodes.begin(), BoundNodes.end(), HasSmallerIndex); + + /*--- Obtaining unique set ---*/ + BoundNodes.resize(std::distance(BoundNodes.begin(), unique(BoundNodes.begin(), BoundNodes.end(), HasEqualIndex))); +} + +void CRadialBasisFunctionInterpolation::SetCtrlNodes(CConfig* config){ + + /*--- Assigning the control nodes based on whether data reduction is applied or not. ---*/ + if(config->GetRBF_DataReduction()){ + + /*--- Control nodes are an empty set ---*/ + ControlNodes = &ReducedControlNodes; + }else{ + + /*--- Control nodes are the boundary nodes ---*/ + ControlNodes = &BoundNodes; + } + + /*--- Obtaining the total number of control nodes. ---*/ + Get_nCtrlNodesGlobal(); + +}; + +void CRadialBasisFunctionInterpolation::ComputeInterpolationMatrix(CGeometry* geometry, const RADIAL_BASIS& type, const su2double radius, su2passivematrix& invInterpMat){ + + /*--- In case of parallel computation, the interpolation coefficients are computed on the master node ---*/ + + if(rank == MASTER_NODE){ + CSymmetricMatrix interpMat; + + /*--- Initialization of the interpolation matrix ---*/ + interpMat.Initialize(nCtrlNodesGlobal); + + /*--- Construction of the interpolation matrix. + Since this matrix is symmetric only upper halve has to be considered ---*/ + + + /*--- Looping over the target nodes ---*/ + for( auto iNode = 0ul; iNode < nCtrlNodesGlobal; iNode++ ){ + + /*--- Looping over the control nodes ---*/ + for ( auto jNode = iNode; jNode < nCtrlNodesGlobal; jNode++){ + + /*--- Distance between nodes ---*/ + auto dist = GeometryToolbox::Distance(nDim, CtrlCoords[iNode*nDim], CtrlCoords[jNode*nDim]); + + /*--- Evaluation of RBF ---*/ + interpMat(iNode, jNode) = SU2_TYPE::GetValue(CRadialBasisFunction::Get_RadialBasisValue(type, radius, dist)); + } + } + + /*--- Obtaining lower halve using symmetry ---*/ + const bool kernelIsSPD = (type == RADIAL_BASIS::WENDLAND_C2) || (type == RADIAL_BASIS::GAUSSIAN) || + (type == RADIAL_BASIS::INV_MULTI_QUADRIC); + + /*--- inverting the interpolation matrix ---*/ + interpMat.Invert(kernelIsSPD); + invInterpMat = interpMat.StealData(); + } +} + +void CRadialBasisFunctionInterpolation::SetDeformation(CGeometry* geometry, CConfig* config){ + + /* --- Initialization of the deformation vector ---*/ + CtrlNodeDeformation.resize(ControlNodes->size()*nDim, 0.0); + + /*--- If requested (no by default) impose the surface deflections in + increments and solve the grid deformation with + successive small deformations. ---*/ + const su2double VarIncrement = 1.0 / ((su2double)config->GetGridDef_Nonlinear_Iter()); + + /*--- Loop over the control nodes ---*/ + for (auto iNode = 0ul; iNode < ControlNodes->size(); iNode++) { + + /*--- Setting nonzero displacement of the moving markers, else setting zero displacement for static markers---*/ + if (config->GetMarker_All_Moving((*ControlNodes)[iNode]->GetMarker())) { + for (auto iDim = 0u; iDim < nDim; iDim++) { + CtrlNodeDeformation[iNode*nDim + iDim] = SU2_TYPE::GetValue(geometry->vertex[(*ControlNodes)[iNode]->GetMarker()][(*ControlNodes)[iNode]->GetVertex()]->GetVarCoord()[iDim] * VarIncrement); + } + } + + else{ + for (auto iDim = 0u; iDim < nDim; iDim++) { + CtrlNodeDeformation[iNode*nDim + iDim] = 0.0; + } + } + } + + /*--- In case of a parallel computation, the deformation of all control nodes is send to the master process ---*/ + #ifdef HAVE_MPI + + /*--- Local number of control nodes ---*/ + unsigned long Local_nControlNodes = ControlNodes->size(); + + /*--- Array containing the local number of control nodes ---*/ + unsigned long Local_nControlNodesArr[size]; + + /*--- gathering local control node coordinate sizes on all processes. ---*/ + SU2_MPI::Allgather(&Local_nControlNodes, 1, MPI_UNSIGNED_LONG, Local_nControlNodesArr, 1, MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); + + /*--- Gathering all deformation vectors on the master node ---*/ + if(rank==MASTER_NODE){ + + /*--- resizing the global deformation vector ---*/ + CtrlNodeDeformation.resize(nCtrlNodesGlobal*nDim); + + /*--- Receiving the local deformation vector from other processes ---*/ + unsigned long start_idx = 0; + for (auto iProc = 0; iProc < size; iProc++) { + if (iProc != MASTER_NODE) { + SU2_MPI::Recv(&CtrlNodeDeformation[0] + start_idx, Local_nControlNodesArr[iProc]*nDim, MPI_DOUBLE, iProc, 0, SU2_MPI::GetComm(), MPI_STATUS_IGNORE); + } + start_idx += Local_nControlNodesArr[iProc]*nDim; + } + + }else{ + + /*--- Sending the local deformation vector to the master node ---*/ + SU2_MPI::Send(CtrlNodeDeformation.data(), Local_nControlNodes*nDim, MPI_DOUBLE, MASTER_NODE, 0, SU2_MPI::GetComm()); + } + #endif +} + +void CRadialBasisFunctionInterpolation::SetInternalNodes(CGeometry* geometry, CConfig* config, vector& internalNodes){ + + /*--- Looping over all nodes and check if part of domain and not on boundary ---*/ + for (auto iNode = 0ul; iNode < geometry->GetnPoint(); iNode++) { + if (!geometry->nodes->GetBoundary(iNode)) { + internalNodes.push_back(iNode); + } + } + + /*--- Adding nodes on markers considered as internal nodes ---*/ + for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); iMarker++){ + + /*--- Check if marker is considered as internal nodes ---*/ + if(config->GetMarker_All_Deform_Mesh_Internal(iMarker)){ + + /*--- Loop over marker vertices ---*/ + for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) { + + /*--- Local node index ---*/ + auto iNode = geometry->vertex[iMarker][iVertex]->GetNode(); + + /*--- if not among the boundary nodes ---*/ + if (find_if (BoundNodes.begin(), BoundNodes.end(), [&](CRadialBasisFunctionNode* i){return i->GetIndex() == iNode;}) == BoundNodes.end()) { + internalNodes.push_back(iNode); + } + } + } + } + + + + /*--- In case of a parallel computation, the nodes on the send/receive markers are included as internal nodes + if they are not already a boundary node with known deformation ---*/ + + #ifdef HAVE_MPI + /*--- Looping over the markers ---*/ + for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); iMarker++) { + + /*--- If send or receive marker ---*/ + if (config->GetMarker_All_SendRecv(iMarker)) { + + /*--- Loop over marker vertices ---*/ + for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) { + + /*--- Local node index ---*/ + auto iNode = geometry->vertex[iMarker][iVertex]->GetNode(); + + /*--- if not among the boundary nodes ---*/ + if (find_if (BoundNodes.begin(), BoundNodes.end(), [&](CRadialBasisFunctionNode* i){return i->GetIndex() == iNode;}) == BoundNodes.end()) { + internalNodes.push_back(iNode); + } + } + } + } + + /*--- sorting of the local indices ---*/ + sort(internalNodes.begin(), internalNodes.end()); + + /*--- Obtaining unique set of internal nodes ---*/ + internalNodes.resize(std::distance(internalNodes.begin(), unique(internalNodes.begin(), internalNodes.end()))); + #endif +} + + +void CRadialBasisFunctionInterpolation::ComputeInterpCoeffs(su2passivematrix& invInterpMat){ + + /*--- resizing the interpolation coefficient vector ---*/ + InterpCoeff.resize(nDim*nCtrlNodesGlobal); + + /*--- Coefficients are found on the master process. + Resulting coefficient is found by summing the multiplications of inverse interpolation matrix entries with deformation ---*/ + if(rank == MASTER_NODE){ + + for(auto iNode = 0ul; iNode < nCtrlNodesGlobal; iNode++){ + for (auto iDim = 0u; iDim < nDim; iDim++){ + InterpCoeff[iNode*nDim + iDim] = 0; + for (auto jNode = 0ul; jNode < nCtrlNodesGlobal; jNode++){ + InterpCoeff[iNode * nDim + iDim] += invInterpMat(iNode,jNode) * CtrlNodeDeformation[jNode*nDim+ iDim]; + } + } + } + } + + /*--- Broadcasting the interpolation coefficients ---*/ + #ifdef HAVE_MPI + SU2_MPI::Bcast(InterpCoeff.data(), InterpCoeff.size(), MPI_DOUBLE, MASTER_NODE, SU2_MPI::GetComm()); + #endif +} + +void CRadialBasisFunctionInterpolation::UpdateGridCoord(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius, const vector& internalNodes){ + + if(rank == MASTER_NODE){ + cout << "updating the grid coordinates" << endl; + } + + /*--- Update of internal node coordinates ---*/ + UpdateInternalCoords(geometry, type, radius, internalNodes); + + /*--- Update of boundary node coordinates ---*/ + UpdateBoundCoords(geometry, config, type, radius); + + /*--- In case of data reduction, perform the correction for nonzero error nodes ---*/ + if(config->GetRBF_DataReduction() && BoundNodes.size() > 0){ + SetCorrection(geometry, config, type, internalNodes); + } +} + +void CRadialBasisFunctionInterpolation::UpdateInternalCoords(CGeometry* geometry, const RADIAL_BASIS& type, const su2double radius, const vector& internalNodes){ + + /*--- Vector for storing the coordinate variation ---*/ + su2double var_coord[nDim]{0.0}; + + /*--- Loop over the internal nodes ---*/ + for(auto iNode = 0ul; iNode < internalNodes.size(); iNode++){ + + /*--- Loop for contribution of each control node ---*/ + for(auto jNode = 0ul; jNode < nCtrlNodesGlobal; jNode++){ + + /*--- Determine distance between considered internal and control node ---*/ + auto dist = GeometryToolbox::Distance(nDim, CtrlCoords[jNode*nDim], geometry->nodes->GetCoord(internalNodes[iNode])); + + /*--- Evaluate RBF based on distance ---*/ + auto rbf = SU2_TYPE::GetValue(CRadialBasisFunction::Get_RadialBasisValue(type, radius, dist)); + + /*--- Add contribution to total coordinate variation ---*/ + for(auto iDim = 0u; iDim < nDim; iDim++){ + var_coord[iDim] += rbf*InterpCoeff[jNode * nDim + iDim]; + } + } + + /*--- Apply the coordinate variation and resetting the var_coord vector to zero ---*/ + for(auto iDim = 0u; iDim < nDim; iDim++){ + geometry->nodes->AddCoord(internalNodes[iNode], iDim, var_coord[iDim]); + var_coord[iDim] = 0; + } + } +} + +void CRadialBasisFunctionInterpolation::UpdateBoundCoords(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius){ + + /*--- Vector for storing the coordinate variation ---*/ + su2double var_coord[nDim]{0.0}; + + /*--- In case of data reduction, the non-control boundary nodes are treated as if they where internal nodes ---*/ + if(config->GetRBF_DataReduction()){ + + /*--- Looping over the non selected boundary nodes ---*/ + for(auto iNode = 0ul; iNode < BoundNodes.size(); iNode++){ + + /*--- Finding contribution of each control node ---*/ + for( auto jNode = 0ul; jNode < nCtrlNodesGlobal; jNode++){ + + /*--- Distance of non-selected boundary node to control node ---*/ + auto dist = GeometryToolbox::Distance(nDim, CtrlCoords[jNode*nDim], geometry->nodes->GetCoord(BoundNodes[iNode]->GetIndex())); + + /*--- Evaluation of the radial basis function based on the distance ---*/ + auto rbf = SU2_TYPE::GetValue(CRadialBasisFunction::Get_RadialBasisValue(type, radius, dist)); + + /*--- Computing and add the resulting coordinate variation ---*/ + for(auto iDim = 0u; iDim < nDim; iDim++){ + var_coord[iDim] += rbf*InterpCoeff[jNode * nDim + iDim]; + } + } + + /*--- Applying the coordinate variation and resetting the var_coord vector*/ + for(auto iDim = 0u; iDim < nDim; iDim++){ + geometry->nodes->AddCoord(BoundNodes[iNode]->GetIndex(), iDim, var_coord[iDim]); + var_coord[iDim] = 0; + } + + } + } + + /*--- Applying the surface deformation, which are stored in the deformation vector ---*/ + for(auto jNode = 0ul; jNode < ControlNodes->size(); jNode++){ + if(config->GetMarker_All_Moving((*ControlNodes)[jNode]->GetMarker())){ + for(auto iDim = 0u; iDim < nDim; iDim++){ + geometry->nodes->AddCoord((*ControlNodes)[jNode]->GetIndex(), iDim, CtrlNodeDeformation[jNode*nDim + iDim]); + } + } + } +} + + +void CRadialBasisFunctionInterpolation::GetInitMaxErrorNode(CGeometry* geometry, CConfig* config, unsigned long& maxErrorNodeLocal, su2double& maxErrorLocal){ + + /*--- Set max error to zero ---*/ + maxErrorLocal = 0.0; + + /*--- Loop over the nodes ---*/ + for(auto iNode = 0ul; iNode < BoundNodes.size(); iNode++){ + + /*--- Compute to squared norm of the deformation ---*/ + su2double normSquaredDeformation = GeometryToolbox::SquaredNorm(nDim, geometry->vertex[BoundNodes[iNode]->GetMarker()][BoundNodes[iNode]->GetVertex()]->GetVarCoord()); + + /*--- In case squared norm deformation is larger than the error, update the error ---*/ + if(normSquaredDeformation > maxErrorLocal){ + maxErrorLocal = normSquaredDeformation; + maxErrorNodeLocal = iNode; + } + } + + /*--- Account for the possibility of applying the deformation in multiple steps ---*/ + maxErrorLocal = sqrt(maxErrorLocal) / ((su2double)config->GetGridDef_Nonlinear_Iter()); +} + + +void CRadialBasisFunctionInterpolation::SetCtrlNodeCoords(CGeometry* geometry){ + /*--- The coordinates of all control nodes are made available on all processes ---*/ + + /*--- resizing the matrix containing the global control node coordinates ---*/ + CtrlCoords.resize(nCtrlNodesGlobal*nDim); + + /*--- Array containing the local control node coordinates ---*/ + su2double localCoords[nDim*ControlNodes->size()]; + + /*--- Storing local control node coordinates ---*/ + for(auto iNode = 0ul; iNode < ControlNodes->size(); iNode++){ + auto coord = geometry->nodes->GetCoord((*ControlNodes)[iNode]->GetIndex()); + for ( auto iDim = 0u ; iDim < nDim; iDim++ ){ + localCoords[ iNode * nDim + iDim ] = coord[iDim]; + } + } + + /*--- Gathering local control node coordinate sizes on all processes. ---*/ + int LocalCoordsSizes[size]; + int localCoordsSize = nDim*ControlNodes->size(); + SU2_MPI::Allgather(&localCoordsSize, 1, MPI_INT, LocalCoordsSizes, 1, MPI_INT, SU2_MPI::GetComm()); + + /*--- Array containing the starting indices for the allgatherv operation */ + int disps[SU2_MPI::GetSize()] = {0}; + + for(auto iProc = 1; iProc < SU2_MPI::GetSize(); iProc++){ + disps[iProc] = disps[iProc-1]+LocalCoordsSizes[iProc-1]; + } + + /*--- Distributing global control node coordinates among all processes ---*/ + SU2_MPI::Allgatherv(&localCoords, localCoordsSize, MPI_DOUBLE, CtrlCoords.data(), LocalCoordsSizes, disps, MPI_DOUBLE, SU2_MPI::GetComm()); +}; + + +void CRadialBasisFunctionInterpolation::GetInterpError(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius, unsigned long& maxErrorNodeLocal, su2double& maxErrorLocal){ + /*--- Array containing the local error ---*/ + su2double localError[nDim]; + + /*--- Magnitude of the local maximum error ---*/ + maxErrorLocal = 0.0; + + /*--- Loop over non-selected boundary nodes ---*/ + for(auto iNode = 0ul; iNode < BoundNodes.size(); iNode++){ + + /*--- Compute nodal error ---*/ + GetNodalError(geometry, config, type, radius, iNode, localError); + + /*--- Setting error ---*/ + BoundNodes[iNode]->SetError(localError, nDim); + + /*--- Compute error magnitude and update local maximum error if necessary ---*/ + su2double errorMagnitude = GeometryToolbox::Norm(nDim, localError); + if(errorMagnitude > maxErrorLocal){ + maxErrorLocal = errorMagnitude; + maxErrorNodeLocal = iNode; + } + } +} + +void CRadialBasisFunctionInterpolation::GetNodalError(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const su2double radius, unsigned long iNode, su2double* localError){ + + /*--- If requested (no by default) impose the surface deflections in increments ---*/ + const su2double VarIncrement = 1.0 / ((su2double)config->GetGridDef_Nonlinear_Iter()); + + /*--- If node is part of a moving boundary then the error is defined as the difference + between the found and prescribed displacements. Thus, here the displacement is substracted from the error ---*/ + if(config->GetMarker_All_Moving(BoundNodes[iNode]->GetMarker())){ + auto displacement = geometry->vertex[BoundNodes[iNode]->GetMarker()][BoundNodes[iNode]->GetVertex()]->GetVarCoord(); + + for(auto iDim = 0u; iDim < nDim; iDim++){ + localError[iDim] = -displacement[iDim] * VarIncrement; + } + }else{ + for(auto iDim = 0u; iDim < nDim; iDim++){ + localError[iDim] = 0; + } + } + + /*--- Resulting displacement from the RBF interpolation is added to the error ---*/ + + /*--- Finding contribution of each control node ---*/ + for(auto jNode = 0ul; jNode < nCtrlNodesGlobal; jNode++){ + + /*--- Distance between non-selected boundary node and control node ---*/ + auto dist = GeometryToolbox::Distance(nDim, CtrlCoords[jNode *nDim], geometry->nodes->GetCoord(BoundNodes[iNode]->GetIndex())); + + /*--- Evaluation of Radial Basis Function ---*/ + auto rbf = SU2_TYPE::GetValue(CRadialBasisFunction::Get_RadialBasisValue(type, radius, dist)); + + /*--- Add contribution to error ---*/ + for(auto iDim = 0u; iDim < nDim; iDim++){ + localError[iDim] += rbf*InterpCoeff[jNode*nDim + iDim]; + } + } +} + +void CRadialBasisFunctionInterpolation::SetCorrection(CGeometry* geometry, CConfig* config, const RADIAL_BASIS& type, const vector& internalNodes){ + + /*--- The non-selected control nodes still have a nonzero error once the maximum error falls below the data reduction tolerance. + This error is applied as correction and interpolated into the volumetric mesh for internal nodes that fall within the correction radius. + To evaluate whether an internal node falls within the correction radius an AD tree is constructed of the boundary nodes, + making it possible to determine the distance to the nearest boundary node. ---*/ + + /*--- Construction of the AD tree consisting of the non-selected boundary nodes ---*/ + + /*--- Number of non-selected boundary nodes ---*/ + const unsigned long nVertexBound = BoundNodes.size(); + + /*--- Vector storing the coordinates of the boundary nodes ---*/ + vector Coord_bound(nDim*nVertexBound); + + /*--- Vector storing the IDs of the boundary nodes ---*/ + vector PointIDs(nVertexBound); + + /*--- Correction Radius, equal to maximum error times a prescribed constant ---*/ + const su2double CorrectionRadius = config->GetRBF_DataRedCorrectionFactor()*MaxErrorGlobal; + + /*--- Storing boundary node information ---*/ + unsigned long i = 0; + unsigned long j = 0; + for(auto iVertex = 0ul; iVertex < nVertexBound; iVertex++){ + auto iNode = BoundNodes[iVertex]->GetIndex(); + PointIDs[i++] = iVertex; + for(auto iDim = 0u; iDim < nDim; iDim++){ + Coord_bound[j++] = geometry->nodes->GetCoord(iNode, iDim); + } + } + + /*--- Construction of AD tree ---*/ + CADTPointsOnlyClass BoundADT(nDim, nVertexBound, Coord_bound.data(), PointIDs.data(), true); + + /*--- ID of nearest boundary node ---*/ + unsigned long pointID; + /*--- Distance to nearest boundary node ---*/ + su2double dist; + /*--- rank of nearest boundary node ---*/ + int rankID; + + /*--- Interpolation of the correction to the internal nodes that fall within the correction radius ---*/ + for(auto iNode = 0ul; iNode < internalNodes.size(); iNode++){ + + /*--- Find nearest node ---*/ + BoundADT.DetermineNearestNode(geometry->nodes->GetCoord(internalNodes[iNode]), dist, pointID, rankID); + + /*--- Get error of nearest node ---*/ + auto err = BoundNodes[pointID]->GetError(); + + /*--- evaluate RBF ---*/ + auto rbf = SU2_TYPE::GetValue(CRadialBasisFunction::Get_RadialBasisValue(type, CorrectionRadius, dist)); + + /*--- Apply correction to the internal node ---*/ + for(auto iDim = 0u; iDim < nDim; iDim++){ + geometry->nodes->AddCoord(internalNodes[iNode], iDim, -rbf*err[iDim]); + } + } + + /*--- Applying the correction to the non-selected boundary nodes ---*/ + for(auto iNode = 0ul; iNode < BoundNodes.size(); iNode++){ + auto err = BoundNodes[iNode]->GetError(); + for(auto iDim = 0u; iDim < nDim; iDim++){ + geometry->nodes->AddCoord(BoundNodes[iNode]->GetIndex(), iDim, -err[iDim]); + } + } +} + + +void CRadialBasisFunctionInterpolation::AddControlNode(unsigned long maxErrorNode){ + /*--- Addition of node to the reduced set of control nodes ---*/ + ReducedControlNodes.push_back(move(BoundNodes[maxErrorNode])); + + /*--- Removal of node among the non-selected boundary nodes ---*/ + BoundNodes.erase(BoundNodes.begin()+maxErrorNode); +} + + +void CRadialBasisFunctionInterpolation::Get_nCtrlNodesGlobal(){ + /*--- Determining the global number of control nodes ---*/ + + /*--- Local number of control nodes ---*/ + auto local_nControlNodes = ControlNodes->size(); + + /*--- Summation of local number of control nodes ---*/ + SU2_MPI::Allreduce(&local_nControlNodes, &nCtrlNodesGlobal, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); +} \ No newline at end of file diff --git a/Common/src/grid_movement/CRadialBasisFunctionNode.cpp b/Common/src/grid_movement/CRadialBasisFunctionNode.cpp new file mode 100644 index 000000000000..4a00aae5ae44 --- /dev/null +++ b/Common/src/grid_movement/CRadialBasisFunctionNode.cpp @@ -0,0 +1,39 @@ +/*! + * \file CRadialBasisFunctionNode.cpp + * \brief Class for defining the nodes in the RBF interpolation. + * \author F. van Steen + * \version 8.0.1 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/grid_movement/CRadialBasisFunctionNode.hpp" + +CRadialBasisFunctionNode::CRadialBasisFunctionNode(unsigned long idx_val, unsigned short marker_val, unsigned long vertex_val){ + /*--- local node index ---*/ + idx = idx_val; + + /*--- local marker index ---*/ + marker_idx = marker_val; + + /*--- local vertex index ---*/ + vertex_idx = vertex_val; +}; \ No newline at end of file diff --git a/Common/src/grid_movement/CVolumetricMovement.cpp b/Common/src/grid_movement/CVolumetricMovement.cpp index 5b09a9cdd716..73ed22099a42 100644 --- a/Common/src/grid_movement/CVolumetricMovement.cpp +++ b/Common/src/grid_movement/CVolumetricMovement.cpp @@ -26,63 +26,22 @@ */ #include "../../include/grid_movement/CVolumetricMovement.hpp" -#include "../../include/adt/CADTPointsOnlyClass.hpp" -#include "../../include/toolboxes/geometry_toolbox.hpp" -CVolumetricMovement::CVolumetricMovement() : CGridMovement(), System(LINEAR_SOLVER_MODE::MESH_DEFORM) {} -CVolumetricMovement::CVolumetricMovement(CGeometry* geometry, CConfig* config) - : CGridMovement(), System(LINEAR_SOLVER_MODE::MESH_DEFORM) { +CVolumetricMovement::CVolumetricMovement() : CGridMovement() {} + +CVolumetricMovement::CVolumetricMovement(CGeometry* geometry) : CGridMovement() { size = SU2_MPI::GetSize(); rank = SU2_MPI::GetRank(); - /*--- Initialize the number of spatial dimensions, length of the state - vector (same as spatial dimensions for grid deformation), and grid nodes. ---*/ - nDim = geometry->GetnDim(); - nVar = geometry->GetnDim(); - nPoint = geometry->GetnPoint(); - nPointDomain = geometry->GetnPointDomain(); - - nIterMesh = 0; - - /*--- Initialize matrix, solution, and r.h.s. structures for the linear solver. ---*/ - if (config->GetVolumetric_Movement() || config->GetSmoothGradient()) { - LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); - LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); - StiffMatrix.Initialize(nPoint, nPointDomain, nVar, nVar, false, geometry, config); - } } CVolumetricMovement::~CVolumetricMovement() = default; -void CVolumetricMovement::UpdateGridCoord(CGeometry* geometry, CConfig* config) { - unsigned short iDim; - unsigned long iPoint, total_index; - su2double new_coord; - - /*--- Update the grid coordinates using the solution of the linear system - after grid deformation (LinSysSol contains the x, y, z displacements). ---*/ - - for (iPoint = 0; iPoint < nPoint; iPoint++) - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - new_coord = geometry->nodes->GetCoord(iPoint, iDim) + LinSysSol[total_index]; - if (fabs(new_coord) < EPS * EPS) new_coord = 0.0; - geometry->nodes->SetCoord(iPoint, iDim, new_coord); - } - - /*--- LinSysSol contains the non-transformed displacements in the periodic halo cells. - * Hence we still need a communication of the transformed coordinates, otherwise periodicity - * is not maintained. ---*/ - - geometry->InitiateComms(geometry, config, COORDINATES); - geometry->CompleteComms(geometry, config, COORDINATES); -} - void CVolumetricMovement::UpdateDualGrid(CGeometry* geometry, CConfig* config) { /*--- After moving all nodes, update the dual mesh. Recompute the edges and - dual mesh control volumes in the domain and on the boundaries. ---*/ +dual mesh control volumes in the domain and on the boundaries. ---*/ geometry->SetControlVolume(config, UPDATE); geometry->SetBoundControlVolume(config, UPDATE); @@ -93,7 +52,7 @@ void CVolumetricMovement::UpdateMultiGrid(CGeometry** geometry, CConfig* config) unsigned short iMGfine, iMGlevel, nMGlevel = config->GetnMGLevels(); /*--- Update the multigrid structure after moving the finest grid, - including computing the grid velocities on the coarser levels. ---*/ +including computing the grid velocities on the coarser levels. ---*/ for (iMGlevel = 1; iMGlevel <= nMGlevel; iMGlevel++) { iMGfine = iMGlevel - 1; @@ -104,121 +63,6 @@ void CVolumetricMovement::UpdateMultiGrid(CGeometry** geometry, CConfig* config) } } -void CVolumetricMovement::SetVolume_Deformation(CGeometry* geometry, CConfig* config, bool UpdateGeo, bool Derivative, - bool ForwardProjectionDerivative) { - unsigned long Tot_Iter = 0; - su2double MinVolume, MaxVolume; - - /*--- Retrieve number or iterations, tol, output, etc. from config ---*/ - - auto Screen_Output = config->GetDeform_Output(); - auto Nonlinear_Iter = config->GetGridDef_Nonlinear_Iter(); - - /*--- Disable the screen output if we're running SU2_CFD ---*/ - - if (config->GetKind_SU2() == SU2_COMPONENT::SU2_CFD && !Derivative) Screen_Output = false; - if (config->GetSmoothGradient()) Screen_Output = true; - - /*--- Set the number of nonlinear iterations to 1 if Derivative computation is enabled ---*/ - - if (Derivative) Nonlinear_Iter = 1; - - /*--- Loop over the total number of grid deformation iterations. The surface - deformation can be divided into increments to help with stability. In - particular, the linear elasticity equations hold only for small deformations. ---*/ - - for (auto iNonlinear_Iter = 0ul; iNonlinear_Iter < Nonlinear_Iter; iNonlinear_Iter++) { - /*--- Initialize vector and sparse matrix ---*/ - - LinSysSol.SetValZero(); - LinSysRes.SetValZero(); - StiffMatrix.SetValZero(); - - /*--- Compute the stiffness matrix entries for all nodes/elements in the - mesh. FEA uses a finite element method discretization of the linear - elasticity equations (transfers element stiffnesses to point-to-point). ---*/ - - MinVolume = SetFEAMethodContributions_Elem(geometry, config); - - /*--- Set the boundary and volume displacements (as prescribed by the - design variable perturbations controlling the surface shape) - as a Dirichlet BC. ---*/ - - SetBoundaryDisplacements(geometry, config); - - /*--- Fix the location of any points in the domain, if requested. ---*/ - - SetDomainDisplacements(geometry, config); - - /*--- Set the boundary derivatives (overrides the actual displacements) ---*/ - - if (Derivative) { - SetBoundaryDerivatives(geometry, config, ForwardProjectionDerivative); - } - - /*--- Communicate any prescribed boundary displacements via MPI, - so that all nodes have the same solution and r.h.s. entries - across all partitions. ---*/ - - CSysMatrixComms::Initiate(LinSysSol, geometry, config); - CSysMatrixComms::Complete(LinSysSol, geometry, config); - - CSysMatrixComms::Initiate(LinSysRes, geometry, config); - CSysMatrixComms::Complete(LinSysRes, geometry, config); - - /*--- Definition of the preconditioner matrix vector multiplication, and linear solver ---*/ - - /*--- To keep legacy behavior ---*/ - System.SetToleranceType(LinearToleranceType::RELATIVE); - - /*--- If we want no derivatives or the direct derivatives, we solve the system using the - * normal matrix vector product and preconditioner. For the mesh sensitivities using - * the discrete adjoint method we solve the system using the transposed matrix. ---*/ - if (!Derivative || ((config->GetKind_SU2() == SU2_COMPONENT::SU2_CFD) && Derivative) || - (config->GetSmoothGradient() && ForwardProjectionDerivative)) { - Tot_Iter = System.Solve(StiffMatrix, LinSysRes, LinSysSol, geometry, config); - - } else if (Derivative && (config->GetKind_SU2() == SU2_COMPONENT::SU2_DOT)) { - Tot_Iter = System.Solve_b(StiffMatrix, LinSysRes, LinSysSol, geometry, config); - } - su2double Residual = System.GetResidual(); - - /*--- Update the grid coordinates and cell volumes using the solution - of the linear system (usol contains the x, y, z displacements). ---*/ - - if (!Derivative) { - UpdateGridCoord(geometry, config); - } else { - UpdateGridCoord_Derivatives(geometry, config, ForwardProjectionDerivative); - } - if (UpdateGeo) { - UpdateDualGrid(geometry, config); - } - - if (!Derivative) { - /*--- Check for failed deformation (negative volumes). ---*/ - - ComputeDeforming_Element_Volume(geometry, MinVolume, MaxVolume, Screen_Output); - - /*--- Calculate amount of nonconvex elements ---*/ - - ComputenNonconvexElements(geometry, Screen_Output); - } - - /*--- Set number of iterations in the mesh update. ---*/ - - Set_nIterMesh(Tot_Iter); - - if (rank == MASTER_NODE && Screen_Output) { - cout << "Non-linear iter.: " << iNonlinear_Iter + 1 << "/" << Nonlinear_Iter << ". Linear iter.: " << Tot_Iter - << ". "; - if (nDim == 2) - cout << "Min. area: " << MinVolume << ". Error: " << Residual << "." << endl; - else - cout << "Min. volume: " << MinVolume << ". Error: " << Residual << "." << endl; - } - } -} void CVolumetricMovement::ComputeDeforming_Element_Volume(CGeometry* geometry, su2double& MinVolume, su2double& MaxVolume, bool Screen_Output) { @@ -364,644 +208,6 @@ void CVolumetricMovement::ComputenNonconvexElements(CGeometry* geometry, bool Sc geometry->SetnNonconvexElements(nNonconvexElements); } -void CVolumetricMovement::ComputeSolid_Wall_Distance(CGeometry* geometry, CConfig* config, su2double& MinDistance, - su2double& MaxDistance) const { - unsigned long nVertex_SolidWall, ii, jj, iVertex, iPoint, pointID; - unsigned short iMarker, iDim; - su2double dist, MaxDistance_Local, MinDistance_Local; - int rankID; - - /*--- Initialize min and max distance ---*/ - - MaxDistance = -1E22; - MinDistance = 1E22; - - /*--- Compute the total number of nodes on no-slip boundaries ---*/ - - nVertex_SolidWall = 0; - for (iMarker = 0; iMarker < config->GetnMarker_All(); ++iMarker) { - if (config->GetSolid_Wall(iMarker)) nVertex_SolidWall += geometry->GetnVertex(iMarker); - } - - /*--- Allocate the vectors to hold boundary node coordinates - and its local ID. ---*/ - - vector Coord_bound(nDim * nVertex_SolidWall); - vector PointIDs(nVertex_SolidWall); - - /*--- Retrieve and store the coordinates of the no-slip boundary nodes - and their local point IDs. ---*/ - - ii = 0; - jj = 0; - for (iMarker = 0; iMarker < config->GetnMarker_All(); ++iMarker) { - if (config->GetSolid_Wall(iMarker)) { - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); ++iVertex) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - PointIDs[jj++] = iPoint; - for (iDim = 0; iDim < nDim; ++iDim) Coord_bound[ii++] = geometry->nodes->GetCoord(iPoint, iDim); - } - } - } - - /*--- Build the ADT of the boundary nodes. ---*/ - - CADTPointsOnlyClass WallADT(nDim, nVertex_SolidWall, Coord_bound.data(), PointIDs.data(), true); - - /*--- Loop over all interior mesh nodes and compute the distances to each - of the no-slip boundary nodes. Store the minimum distance to the wall - for each interior mesh node. ---*/ - - if (WallADT.IsEmpty()) { - /*--- No solid wall boundary nodes in the entire mesh. - Set the wall distance to zero for all nodes. ---*/ - - for (iPoint = 0; iPoint < geometry->GetnPoint(); ++iPoint) geometry->nodes->SetWall_Distance(iPoint, 0.0); - } else { - /*--- Solid wall boundary nodes are present. Compute the wall - distance for all nodes. ---*/ - - for (iPoint = 0; iPoint < geometry->GetnPoint(); ++iPoint) { - WallADT.DetermineNearestNode(geometry->nodes->GetCoord(iPoint), dist, pointID, rankID); - geometry->nodes->SetWall_Distance(iPoint, dist); - - MaxDistance = max(MaxDistance, dist); - - /*--- To discard points on the surface we use > EPS ---*/ - - if (sqrt(dist) > EPS) MinDistance = min(MinDistance, dist); - } - - MaxDistance_Local = MaxDistance; - MaxDistance = 0.0; - MinDistance_Local = MinDistance; - MinDistance = 0.0; - -#ifdef HAVE_MPI - SU2_MPI::Allreduce(&MaxDistance_Local, &MaxDistance, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm()); - SU2_MPI::Allreduce(&MinDistance_Local, &MinDistance, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm()); -#else - MaxDistance = MaxDistance_Local; - MinDistance = MinDistance_Local; -#endif - } -} - -su2double CVolumetricMovement::SetFEAMethodContributions_Elem(CGeometry* geometry, CConfig* config) { - unsigned short iVar, iDim, nNodes = 0, iNodes, StiffMatrix_nElem = 0; - unsigned long iElem, PointCorners[8]; - su2double **StiffMatrix_Elem = nullptr, CoordCorners[8][3]; - su2double MinVolume = 0.0, MaxVolume = 0.0, MinDistance = 0.0, MaxDistance = 0.0, ElemVolume = 0.0, - ElemDistance = 0.0; - - bool Screen_Output = config->GetDeform_Output(); - - /*--- Allocate maximum size (quadrilateral and hexahedron) ---*/ - - if (nDim == 2) - StiffMatrix_nElem = 8; - else - StiffMatrix_nElem = 24; - - StiffMatrix_Elem = new su2double*[StiffMatrix_nElem]; - for (iVar = 0; iVar < StiffMatrix_nElem; iVar++) StiffMatrix_Elem[iVar] = new su2double[StiffMatrix_nElem]; - - /*--- Compute min volume in the entire mesh. ---*/ - - ComputeDeforming_Element_Volume(geometry, MinVolume, MaxVolume, Screen_Output); - if (rank == MASTER_NODE && Screen_Output) - cout << "Min. volume: " << MinVolume << ", max. volume: " << MaxVolume << "." << endl; - - /*--- Compute the distance to the nearest surface if needed - as part of the stiffness calculation.. ---*/ - - if ((config->GetDeform_Stiffness_Type() == SOLID_WALL_DISTANCE) || (config->GetDeform_Limit() < 1E6)) { - ComputeSolid_Wall_Distance(geometry, config, MinDistance, MaxDistance); - if (rank == MASTER_NODE && Screen_Output) - cout << "Min. distance: " << MinDistance << ", max. distance: " << MaxDistance << "." << endl; - } - - /*--- Compute contributions from each element by forming the stiffness matrix (FEA) ---*/ - - for (iElem = 0; iElem < geometry->GetnElem(); iElem++) { - if (geometry->elem[iElem]->GetVTK_Type() == TRIANGLE) nNodes = 3; - if (geometry->elem[iElem]->GetVTK_Type() == QUADRILATERAL) nNodes = 4; - if (geometry->elem[iElem]->GetVTK_Type() == TETRAHEDRON) nNodes = 4; - if (geometry->elem[iElem]->GetVTK_Type() == PYRAMID) nNodes = 5; - if (geometry->elem[iElem]->GetVTK_Type() == PRISM) nNodes = 6; - if (geometry->elem[iElem]->GetVTK_Type() == HEXAHEDRON) nNodes = 8; - - for (iNodes = 0; iNodes < nNodes; iNodes++) { - PointCorners[iNodes] = geometry->elem[iElem]->GetNode(iNodes); - for (iDim = 0; iDim < nDim; iDim++) { - CoordCorners[iNodes][iDim] = geometry->nodes->GetCoord(PointCorners[iNodes], iDim); - } - } - - /*--- Extract Element volume and distance to compute the stiffness ---*/ - - ElemVolume = geometry->elem[iElem]->GetVolume(); - - if ((config->GetDeform_Stiffness_Type() == SOLID_WALL_DISTANCE)) { - ElemDistance = 0.0; - for (iNodes = 0; iNodes < nNodes; iNodes++) - ElemDistance += geometry->nodes->GetWall_Distance(PointCorners[iNodes]); - ElemDistance = ElemDistance / (su2double)nNodes; - } - - if (nDim == 2) - SetFEA_StiffMatrix2D(geometry, config, StiffMatrix_Elem, PointCorners, CoordCorners, nNodes, ElemVolume, - ElemDistance); - if (nDim == 3) - SetFEA_StiffMatrix3D(geometry, config, StiffMatrix_Elem, PointCorners, CoordCorners, nNodes, ElemVolume, - ElemDistance); - - AddFEA_StiffMatrix(geometry, StiffMatrix_Elem, PointCorners, nNodes); - } - - /*--- Deallocate memory and exit ---*/ - - for (iVar = 0; iVar < StiffMatrix_nElem; iVar++) delete[] StiffMatrix_Elem[iVar]; - delete[] StiffMatrix_Elem; - - return MinVolume; -} - -su2double CVolumetricMovement::ShapeFunc_Triangle(su2double Xi, su2double Eta, su2double CoordCorners[8][3], - su2double DShapeFunction[8][4]) { - int i, j, k; - su2double c0, c1, xsj; - su2double xs[3][3], ad[3][3]; - - /*--- Shape functions ---*/ - - DShapeFunction[0][3] = Xi; - DShapeFunction[1][3] = Eta; - DShapeFunction[2][3] = 1 - Xi - Eta; - - /*--- dN/d xi, dN/d eta ---*/ - - DShapeFunction[0][0] = 1.0; - DShapeFunction[0][1] = 0.0; - DShapeFunction[1][0] = 0.0; - DShapeFunction[1][1] = 1.0; - DShapeFunction[2][0] = -1.0; - DShapeFunction[2][1] = -1.0; - - /*--- Jacobian transformation ---*/ - - for (i = 0; i < 2; i++) { - for (j = 0; j < 2; j++) { - xs[i][j] = 0.0; - for (k = 0; k < 3; k++) { - xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; - } - } - } - - /*--- Adjoint to Jacobian ---*/ - - ad[0][0] = xs[1][1]; - ad[0][1] = -xs[0][1]; - ad[1][0] = -xs[1][0]; - ad[1][1] = xs[0][0]; - - /*--- Determinant of Jacobian ---*/ - - xsj = ad[0][0] * ad[1][1] - ad[0][1] * ad[1][0]; - - /*--- Jacobian inverse ---*/ - - for (i = 0; i < 2; i++) { - for (j = 0; j < 2; j++) { - xs[i][j] = ad[i][j] / xsj; - } - } - - /*--- Derivatives with repect to global coordinates ---*/ - - for (k = 0; k < 3; k++) { - c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1]; // dN/dx - c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1]; // dN/dy - DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi - DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta - } - - return xsj; -} - -su2double CVolumetricMovement::ShapeFunc_Quadrilateral(su2double Xi, su2double Eta, su2double CoordCorners[8][3], - su2double DShapeFunction[8][4]) { - int i, j, k; - su2double c0, c1, xsj; - su2double xs[3][3], ad[3][3]; - - /*--- Shape functions ---*/ - - DShapeFunction[0][3] = 0.25 * (1.0 - Xi) * (1.0 - Eta); - DShapeFunction[1][3] = 0.25 * (1.0 + Xi) * (1.0 - Eta); - DShapeFunction[2][3] = 0.25 * (1.0 + Xi) * (1.0 + Eta); - DShapeFunction[3][3] = 0.25 * (1.0 - Xi) * (1.0 + Eta); - - /*--- dN/d xi, dN/d eta ---*/ - - DShapeFunction[0][0] = -0.25 * (1.0 - Eta); - DShapeFunction[0][1] = -0.25 * (1.0 - Xi); - DShapeFunction[1][0] = 0.25 * (1.0 - Eta); - DShapeFunction[1][1] = -0.25 * (1.0 + Xi); - DShapeFunction[2][0] = 0.25 * (1.0 + Eta); - DShapeFunction[2][1] = 0.25 * (1.0 + Xi); - DShapeFunction[3][0] = -0.25 * (1.0 + Eta); - DShapeFunction[3][1] = 0.25 * (1.0 - Xi); - - /*--- Jacobian transformation ---*/ - - for (i = 0; i < 2; i++) { - for (j = 0; j < 2; j++) { - xs[i][j] = 0.0; - for (k = 0; k < 4; k++) { - xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; - } - } - } - - /*--- Adjoint to Jacobian ---*/ - - ad[0][0] = xs[1][1]; - ad[0][1] = -xs[0][1]; - ad[1][0] = -xs[1][0]; - ad[1][1] = xs[0][0]; - - /*--- Determinant of Jacobian ---*/ - - xsj = ad[0][0] * ad[1][1] - ad[0][1] * ad[1][0]; - - /*--- Jacobian inverse ---*/ - - for (i = 0; i < 2; i++) { - for (j = 0; j < 2; j++) { - xs[i][j] = ad[i][j] / xsj; - } - } - - /*--- Derivatives with repect to global coordinates ---*/ - - for (k = 0; k < 4; k++) { - c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1]; // dN/dx - c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1]; // dN/dy - DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi - DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta - } - - return xsj; -} - -su2double CVolumetricMovement::ShapeFunc_Tetra(su2double Xi, su2double Eta, su2double Zeta, - su2double CoordCorners[8][3], su2double DShapeFunction[8][4]) { - int i, j, k; - su2double c0, c1, c2, xsj; - su2double xs[3][3], ad[3][3]; - - /*--- Shape functions ---*/ - - DShapeFunction[0][3] = Xi; - DShapeFunction[1][3] = Zeta; - DShapeFunction[2][3] = 1.0 - Xi - Eta - Zeta; - DShapeFunction[3][3] = Eta; - - /*--- dN/d xi, dN/d eta, dN/d zeta ---*/ - - DShapeFunction[0][0] = 1.0; - DShapeFunction[0][1] = 0.0; - DShapeFunction[0][2] = 0.0; - DShapeFunction[1][0] = 0.0; - DShapeFunction[1][1] = 0.0; - DShapeFunction[1][2] = 1.0; - DShapeFunction[2][0] = -1.0; - DShapeFunction[2][1] = -1.0; - DShapeFunction[2][2] = -1.0; - DShapeFunction[3][0] = 0.0; - DShapeFunction[3][1] = 1.0; - DShapeFunction[3][2] = 0.0; - - /*--- Jacobian transformation ---*/ - - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - xs[i][j] = 0.0; - for (k = 0; k < 4; k++) { - xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; - } - } - } - - /*--- Adjoint to Jacobian ---*/ - - ad[0][0] = xs[1][1] * xs[2][2] - xs[1][2] * xs[2][1]; - ad[0][1] = xs[0][2] * xs[2][1] - xs[0][1] * xs[2][2]; - ad[0][2] = xs[0][1] * xs[1][2] - xs[0][2] * xs[1][1]; - ad[1][0] = xs[1][2] * xs[2][0] - xs[1][0] * xs[2][2]; - ad[1][1] = xs[0][0] * xs[2][2] - xs[0][2] * xs[2][0]; - ad[1][2] = xs[0][2] * xs[1][0] - xs[0][0] * xs[1][2]; - ad[2][0] = xs[1][0] * xs[2][1] - xs[1][1] * xs[2][0]; - ad[2][1] = xs[0][1] * xs[2][0] - xs[0][0] * xs[2][1]; - ad[2][2] = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; - - /*--- Determinant of Jacobian ---*/ - - xsj = xs[0][0] * ad[0][0] + xs[0][1] * ad[1][0] + xs[0][2] * ad[2][0]; - - /*--- Jacobian inverse ---*/ - - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - xs[i][j] = ad[i][j] / xsj; - } - } - - /*--- Derivatives with repect to global coordinates ---*/ - - for (k = 0; k < 4; k++) { - c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1] + xs[0][2] * DShapeFunction[k][2]; // dN/dx - c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1] + xs[1][2] * DShapeFunction[k][2]; // dN/dy - c2 = xs[2][0] * DShapeFunction[k][0] + xs[2][1] * DShapeFunction[k][1] + xs[2][2] * DShapeFunction[k][2]; // dN/dz - DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi - DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta - DShapeFunction[k][2] = c2; // store dN/dz instead of dN/d zeta - } - - return xsj; -} - -su2double CVolumetricMovement::ShapeFunc_Pyram(su2double Xi, su2double Eta, su2double Zeta, - su2double CoordCorners[8][3], su2double DShapeFunction[8][4]) { - int i, j, k; - su2double c0, c1, c2, xsj; - su2double xs[3][3], ad[3][3]; - - /*--- Shape functions ---*/ - - DShapeFunction[0][3] = 0.25 * (-Xi + Eta + Zeta - 1.0) * (-Xi - Eta + Zeta - 1.0) / (1.0 - Zeta); - DShapeFunction[1][3] = 0.25 * (-Xi - Eta + Zeta - 1.0) * (Xi - Eta + Zeta - 1.0) / (1.0 - Zeta); - DShapeFunction[2][3] = 0.25 * (Xi + Eta + Zeta - 1.0) * (Xi - Eta + Zeta - 1.0) / (1.0 - Zeta); - DShapeFunction[3][3] = 0.25 * (Xi + Eta + Zeta - 1.0) * (-Xi + Eta + Zeta - 1.0) / (1.0 - Zeta); - DShapeFunction[4][3] = Zeta; - - /*--- dN/d xi ---*/ - - DShapeFunction[0][0] = 0.5 * (Zeta - Xi - 1.0) / (Zeta - 1.0); - DShapeFunction[1][0] = 0.5 * Xi / (Zeta - 1.0); - DShapeFunction[2][0] = 0.5 * (1.0 - Zeta - Xi) / (Zeta - 1.0); - DShapeFunction[3][0] = DShapeFunction[1][0]; - DShapeFunction[4][0] = 0.0; - - /*--- dN/d eta ---*/ - - DShapeFunction[0][1] = 0.5 * Eta / (Zeta - 1.0); - DShapeFunction[1][1] = 0.5 * (Zeta - Eta - 1.0) / (Zeta - 1.0); - DShapeFunction[2][1] = DShapeFunction[0][1]; - DShapeFunction[3][1] = 0.5 * (1.0 - Zeta - Eta) / (Zeta - 1.0); - DShapeFunction[4][1] = 0.0; - - /*--- dN/d zeta ---*/ - - DShapeFunction[0][2] = 0.25 * (-1.0 + 2.0 * Zeta - Zeta * Zeta - Eta * Eta + Xi * Xi) / ((1.0 - Zeta) * (1.0 - Zeta)); - DShapeFunction[1][2] = 0.25 * (-1.0 + 2.0 * Zeta - Zeta * Zeta + Eta * Eta - Xi * Xi) / ((1.0 - Zeta) * (1.0 - Zeta)); - DShapeFunction[2][2] = DShapeFunction[0][2]; - DShapeFunction[3][2] = DShapeFunction[1][2]; - DShapeFunction[4][2] = 1.0; - - /*--- Jacobian transformation ---*/ - - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - xs[i][j] = 0.0; - for (k = 0; k < 5; k++) { - xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; - } - } - } - - /*--- Adjoint to Jacobian ---*/ - - ad[0][0] = xs[1][1] * xs[2][2] - xs[1][2] * xs[2][1]; - ad[0][1] = xs[0][2] * xs[2][1] - xs[0][1] * xs[2][2]; - ad[0][2] = xs[0][1] * xs[1][2] - xs[0][2] * xs[1][1]; - ad[1][0] = xs[1][2] * xs[2][0] - xs[1][0] * xs[2][2]; - ad[1][1] = xs[0][0] * xs[2][2] - xs[0][2] * xs[2][0]; - ad[1][2] = xs[0][2] * xs[1][0] - xs[0][0] * xs[1][2]; - ad[2][0] = xs[1][0] * xs[2][1] - xs[1][1] * xs[2][0]; - ad[2][1] = xs[0][1] * xs[2][0] - xs[0][0] * xs[2][1]; - ad[2][2] = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; - - /*--- Determinant of Jacobian ---*/ - - xsj = xs[0][0] * ad[0][0] + xs[0][1] * ad[1][0] + xs[0][2] * ad[2][0]; - - /*--- Jacobian inverse ---*/ - - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - xs[i][j] = ad[i][j] / xsj; - } - } - - /*--- Derivatives with repect to global coordinates ---*/ - - for (k = 0; k < 5; k++) { - c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1] + xs[0][2] * DShapeFunction[k][2]; // dN/dx - c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1] + xs[1][2] * DShapeFunction[k][2]; // dN/dy - c2 = xs[2][0] * DShapeFunction[k][0] + xs[2][1] * DShapeFunction[k][1] + xs[2][2] * DShapeFunction[k][2]; // dN/dz - DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi - DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta - DShapeFunction[k][2] = c2; // store dN/dz instead of dN/d zeta - } - - return xsj; -} - -su2double CVolumetricMovement::ShapeFunc_Prism(su2double Xi, su2double Eta, su2double Zeta, - su2double CoordCorners[8][3], su2double DShapeFunction[8][4]) { - int i, j, k; - su2double c0, c1, c2, xsj; - su2double xs[3][3], ad[3][3]; - - /*--- Shape functions ---*/ - - DShapeFunction[0][3] = 0.5 * Eta * (1.0 - Xi); - DShapeFunction[1][3] = 0.5 * Zeta * (1.0 - Xi); - DShapeFunction[2][3] = 0.5 * (1.0 - Eta - Zeta) * (1.0 - Xi); - DShapeFunction[3][3] = 0.5 * Eta * (Xi + 1.0); - DShapeFunction[4][3] = 0.5 * Zeta * (Xi + 1.0); - DShapeFunction[5][3] = 0.5 * (1.0 - Eta - Zeta) * (Xi + 1.0); - - /*--- dN/d Xi, dN/d Eta, dN/d Zeta ---*/ - - DShapeFunction[0][0] = -0.5 * Eta; - DShapeFunction[0][1] = 0.5 * (1.0 - Xi); - DShapeFunction[0][2] = 0.0; - DShapeFunction[1][0] = -0.5 * Zeta; - DShapeFunction[1][1] = 0.0; - DShapeFunction[1][2] = 0.5 * (1.0 - Xi); - DShapeFunction[2][0] = -0.5 * (1.0 - Eta - Zeta); - DShapeFunction[2][1] = -0.5 * (1.0 - Xi); - DShapeFunction[2][2] = -0.5 * (1.0 - Xi); - DShapeFunction[3][0] = 0.5 * Eta; - DShapeFunction[3][1] = 0.5 * (Xi + 1.0); - DShapeFunction[3][2] = 0.0; - DShapeFunction[4][0] = 0.5 * Zeta; - DShapeFunction[4][1] = 0.0; - DShapeFunction[4][2] = 0.5 * (Xi + 1.0); - DShapeFunction[5][0] = 0.5 * (1.0 - Eta - Zeta); - DShapeFunction[5][1] = -0.5 * (Xi + 1.0); - DShapeFunction[5][2] = -0.5 * (Xi + 1.0); - - /*--- Jacobian transformation ---*/ - - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - xs[i][j] = 0.0; - for (k = 0; k < 6; k++) { - xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; - } - } - } - - /*--- Adjoint to Jacobian ---*/ - - ad[0][0] = xs[1][1] * xs[2][2] - xs[1][2] * xs[2][1]; - ad[0][1] = xs[0][2] * xs[2][1] - xs[0][1] * xs[2][2]; - ad[0][2] = xs[0][1] * xs[1][2] - xs[0][2] * xs[1][1]; - ad[1][0] = xs[1][2] * xs[2][0] - xs[1][0] * xs[2][2]; - ad[1][1] = xs[0][0] * xs[2][2] - xs[0][2] * xs[2][0]; - ad[1][2] = xs[0][2] * xs[1][0] - xs[0][0] * xs[1][2]; - ad[2][0] = xs[1][0] * xs[2][1] - xs[1][1] * xs[2][0]; - ad[2][1] = xs[0][1] * xs[2][0] - xs[0][0] * xs[2][1]; - ad[2][2] = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; - - /*--- Determinant of Jacobian ---*/ - - xsj = xs[0][0] * ad[0][0] + xs[0][1] * ad[1][0] + xs[0][2] * ad[2][0]; - - /*--- Jacobian inverse ---*/ - - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - xs[i][j] = ad[i][j] / xsj; - } - } - - /*--- Derivatives with repect to global coordinates ---*/ - - for (k = 0; k < 6; k++) { - c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1] + xs[0][2] * DShapeFunction[k][2]; // dN/dx - c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1] + xs[1][2] * DShapeFunction[k][2]; // dN/dy - c2 = xs[2][0] * DShapeFunction[k][0] + xs[2][1] * DShapeFunction[k][1] + xs[2][2] * DShapeFunction[k][2]; // dN/dz - DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi - DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta - DShapeFunction[k][2] = c2; // store dN/dz instead of dN/d zeta - } - - return xsj; -} - -su2double CVolumetricMovement::ShapeFunc_Hexa(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3], - su2double DShapeFunction[8][4]) { - int i, j, k; - su2double c0, c1, c2, xsj; - su2double xs[3][3], ad[3][3]; - - /*--- Shape functions ---*/ - - DShapeFunction[0][3] = 0.125 * (1.0 - Xi) * (1.0 - Eta) * (1.0 - Zeta); - DShapeFunction[1][3] = 0.125 * (1.0 + Xi) * (1.0 - Eta) * (1.0 - Zeta); - DShapeFunction[2][3] = 0.125 * (1.0 + Xi) * (1.0 + Eta) * (1.0 - Zeta); - DShapeFunction[3][3] = 0.125 * (1.0 - Xi) * (1.0 + Eta) * (1.0 - Zeta); - DShapeFunction[4][3] = 0.125 * (1.0 - Xi) * (1.0 - Eta) * (1.0 + Zeta); - DShapeFunction[5][3] = 0.125 * (1.0 + Xi) * (1.0 - Eta) * (1.0 + Zeta); - DShapeFunction[6][3] = 0.125 * (1.0 + Xi) * (1.0 + Eta) * (1.0 + Zeta); - DShapeFunction[7][3] = 0.125 * (1.0 - Xi) * (1.0 + Eta) * (1.0 + Zeta); - - /*--- dN/d xi ---*/ - - DShapeFunction[0][0] = -0.125 * (1.0 - Eta) * (1.0 - Zeta); - DShapeFunction[1][0] = 0.125 * (1.0 - Eta) * (1.0 - Zeta); - DShapeFunction[2][0] = 0.125 * (1.0 + Eta) * (1.0 - Zeta); - DShapeFunction[3][0] = -0.125 * (1.0 + Eta) * (1.0 - Zeta); - DShapeFunction[4][0] = -0.125 * (1.0 - Eta) * (1.0 + Zeta); - DShapeFunction[5][0] = 0.125 * (1.0 - Eta) * (1.0 + Zeta); - DShapeFunction[6][0] = 0.125 * (1.0 + Eta) * (1.0 + Zeta); - DShapeFunction[7][0] = -0.125 * (1.0 + Eta) * (1.0 + Zeta); - - /*--- dN/d eta ---*/ - - DShapeFunction[0][1] = -0.125 * (1.0 - Xi) * (1.0 - Zeta); - DShapeFunction[1][1] = -0.125 * (1.0 + Xi) * (1.0 - Zeta); - DShapeFunction[2][1] = 0.125 * (1.0 + Xi) * (1.0 - Zeta); - DShapeFunction[3][1] = 0.125 * (1.0 - Xi) * (1.0 - Zeta); - DShapeFunction[4][1] = -0.125 * (1.0 - Xi) * (1.0 + Zeta); - DShapeFunction[5][1] = -0.125 * (1.0 + Xi) * (1.0 + Zeta); - DShapeFunction[6][1] = 0.125 * (1.0 + Xi) * (1.0 + Zeta); - DShapeFunction[7][1] = 0.125 * (1.0 - Xi) * (1.0 + Zeta); - - /*--- dN/d zeta ---*/ - - DShapeFunction[0][2] = -0.125 * (1.0 - Xi) * (1.0 - Eta); - DShapeFunction[1][2] = -0.125 * (1.0 + Xi) * (1.0 - Eta); - DShapeFunction[2][2] = -0.125 * (1.0 + Xi) * (1.0 + Eta); - DShapeFunction[3][2] = -0.125 * (1.0 - Xi) * (1.0 + Eta); - DShapeFunction[4][2] = 0.125 * (1.0 - Xi) * (1.0 - Eta); - DShapeFunction[5][2] = 0.125 * (1.0 + Xi) * (1.0 - Eta); - DShapeFunction[6][2] = 0.125 * (1.0 + Xi) * (1.0 + Eta); - DShapeFunction[7][2] = 0.125 * (1.0 - Xi) * (1.0 + Eta); - - /*--- Jacobian transformation ---*/ - - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - xs[i][j] = 0.0; - for (k = 0; k < 8; k++) { - xs[i][j] = xs[i][j] + CoordCorners[k][j] * DShapeFunction[k][i]; - } - } - } - - /*--- Adjoint to Jacobian ---*/ - - ad[0][0] = xs[1][1] * xs[2][2] - xs[1][2] * xs[2][1]; - ad[0][1] = xs[0][2] * xs[2][1] - xs[0][1] * xs[2][2]; - ad[0][2] = xs[0][1] * xs[1][2] - xs[0][2] * xs[1][1]; - ad[1][0] = xs[1][2] * xs[2][0] - xs[1][0] * xs[2][2]; - ad[1][1] = xs[0][0] * xs[2][2] - xs[0][2] * xs[2][0]; - ad[1][2] = xs[0][2] * xs[1][0] - xs[0][0] * xs[1][2]; - ad[2][0] = xs[1][0] * xs[2][1] - xs[1][1] * xs[2][0]; - ad[2][1] = xs[0][1] * xs[2][0] - xs[0][0] * xs[2][1]; - ad[2][2] = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; - - /*--- Determinant of Jacobian ---*/ - - xsj = xs[0][0] * ad[0][0] + xs[0][1] * ad[1][0] + xs[0][2] * ad[2][0]; - - /*--- Jacobian inverse ---*/ - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - xs[i][j] = ad[i][j] / xsj; - } - } - - /*--- Derivatives with repect to global coordinates ---*/ - - for (k = 0; k < 8; k++) { - c0 = xs[0][0] * DShapeFunction[k][0] + xs[0][1] * DShapeFunction[k][1] + xs[0][2] * DShapeFunction[k][2]; // dN/dx - c1 = xs[1][0] * DShapeFunction[k][0] + xs[1][1] * DShapeFunction[k][1] + xs[1][2] * DShapeFunction[k][2]; // dN/dy - c2 = xs[2][0] * DShapeFunction[k][0] + xs[2][1] * DShapeFunction[k][1] + xs[2][2] * DShapeFunction[k][2]; // dN/dz - DShapeFunction[k][0] = c0; // store dN/dx instead of dN/d xi - DShapeFunction[k][1] = c1; // store dN/dy instead of dN/d eta - DShapeFunction[k][2] = c2; // store dN/dz instead of dN/d zeta - } - - return xsj; -} su2double CVolumetricMovement::GetTriangle_Area(su2double CoordCorners[8][3]) const { unsigned short iDim; @@ -1275,651 +481,8 @@ su2double CVolumetricMovement::GetHexa_Volume(su2double CoordCorners[8][3]) cons return Volume; } -void CVolumetricMovement::SetFEA_StiffMatrix2D(CGeometry* geometry, CConfig* config, su2double** StiffMatrix_Elem, - unsigned long PointCorners[8], su2double CoordCorners[8][3], - unsigned short nNodes, su2double ElemVolume, su2double ElemDistance) { - su2double B_Matrix[3][8], D_Matrix[3][3], Aux_Matrix[8][3]; - su2double Xi = 0.0, Eta = 0.0, Det = 0.0, E = 1 / EPS, Lambda = 0.0, Mu = 0.0, Nu = 0.0; - unsigned short iNode, iVar, jVar, kVar, iGauss, nGauss = 0; - su2double DShapeFunction[8][4] = {{0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}}; - su2double Location[4][3], Weight[4]; - unsigned short nVar = geometry->GetnDim(); - - for (iVar = 0; iVar < nNodes * nVar; iVar++) { - for (jVar = 0; jVar < nNodes * nVar; jVar++) { - StiffMatrix_Elem[iVar][jVar] = 0.0; - } - } - - /*--- Integration formulae from "Shape functions and points of - integration of the Résumé" by Josselin DELMAS (2013) ---*/ - - /*--- Triangle. Nodes of numerical integration at 1 point (order 1). ---*/ - - if (nNodes == 3) { - nGauss = 1; - Location[0][0] = 0.333333333333333; - Location[0][1] = 0.333333333333333; - Weight[0] = 0.5; - } - - /*--- Quadrilateral. Nodes of numerical integration at 4 points (order 2). ---*/ - - if (nNodes == 4) { - nGauss = 4; - Location[0][0] = -0.577350269189626; - Location[0][1] = -0.577350269189626; - Weight[0] = 1.0; - Location[1][0] = 0.577350269189626; - Location[1][1] = -0.577350269189626; - Weight[1] = 1.0; - Location[2][0] = 0.577350269189626; - Location[2][1] = 0.577350269189626; - Weight[2] = 1.0; - Location[3][0] = -0.577350269189626; - Location[3][1] = 0.577350269189626; - Weight[3] = 1.0; - } - - for (iGauss = 0; iGauss < nGauss; iGauss++) { - Xi = Location[iGauss][0]; - Eta = Location[iGauss][1]; - - if (nNodes == 3) Det = ShapeFunc_Triangle(Xi, Eta, CoordCorners, DShapeFunction); - if (nNodes == 4) Det = ShapeFunc_Quadrilateral(Xi, Eta, CoordCorners, DShapeFunction); - - /*--- Compute the B Matrix ---*/ - - for (iVar = 0; iVar < 3; iVar++) - for (jVar = 0; jVar < nNodes * nVar; jVar++) B_Matrix[iVar][jVar] = 0.0; - - for (iNode = 0; iNode < nNodes; iNode++) { - B_Matrix[0][0 + iNode * nVar] = DShapeFunction[iNode][0]; - B_Matrix[1][1 + iNode * nVar] = DShapeFunction[iNode][1]; - - B_Matrix[2][0 + iNode * nVar] = DShapeFunction[iNode][1]; - B_Matrix[2][1 + iNode * nVar] = DShapeFunction[iNode][0]; - } - - /*--- Impose a type of stiffness for each element ---*/ - - switch (config->GetDeform_Stiffness_Type()) { - case INVERSE_VOLUME: - E = 1.0 / ElemVolume; - break; - case SOLID_WALL_DISTANCE: - E = 1.0 / ElemDistance; - break; - case CONSTANT_STIFFNESS: - E = 1.0 / EPS; - break; - } - - Nu = config->GetDeform_Coeff(); - Mu = E / (2.0 * (1.0 + Nu)); - Lambda = Nu * E / ((1.0 + Nu) * (1.0 - 2.0 * Nu)); - - /*--- Compute the D Matrix (for plane strain and 3-D)---*/ - - D_Matrix[0][0] = Lambda + 2.0 * Mu; - D_Matrix[0][1] = Lambda; - D_Matrix[0][2] = 0.0; - D_Matrix[1][0] = Lambda; - D_Matrix[1][1] = Lambda + 2.0 * Mu; - D_Matrix[1][2] = 0.0; - D_Matrix[2][0] = 0.0; - D_Matrix[2][1] = 0.0; - D_Matrix[2][2] = Mu; - - /*--- Compute the BT.D Matrix ---*/ - - for (iVar = 0; iVar < nNodes * nVar; iVar++) { - for (jVar = 0; jVar < 3; jVar++) { - Aux_Matrix[iVar][jVar] = 0.0; - for (kVar = 0; kVar < 3; kVar++) Aux_Matrix[iVar][jVar] += B_Matrix[kVar][iVar] * D_Matrix[kVar][jVar]; - } - } - - /*--- Compute the BT.D.B Matrix (stiffness matrix), and add to the original - matrix using Gauss integration ---*/ - - for (iVar = 0; iVar < nNodes * nVar; iVar++) { - for (jVar = 0; jVar < nNodes * nVar; jVar++) { - for (kVar = 0; kVar < 3; kVar++) { - StiffMatrix_Elem[iVar][jVar] += Weight[iGauss] * Aux_Matrix[iVar][kVar] * B_Matrix[kVar][jVar] * fabs(Det); - } - } - } - } -} - -void CVolumetricMovement::SetFEA_StiffMatrix3D(CGeometry* geometry, CConfig* config, su2double** StiffMatrix_Elem, - unsigned long PointCorners[8], su2double CoordCorners[8][3], - unsigned short nNodes, su2double ElemVolume, su2double ElemDistance) { - su2double B_Matrix[6][24], - D_Matrix[6][6] = {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}, - Aux_Matrix[24][6]; - su2double Xi = 0.0, Eta = 0.0, Zeta = 0.0, Det = 0.0, Mu = 0.0, E = 0.0, Lambda = 0.0, Nu = 0.0; - unsigned short iNode, iVar, jVar, kVar, iGauss, nGauss = 0; - su2double DShapeFunction[8][4] = {{0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0}}; - su2double Location[8][3], Weight[8]; - unsigned short nVar = geometry->GetnDim(); - - for (iVar = 0; iVar < nNodes * nVar; iVar++) { - for (jVar = 0; jVar < nNodes * nVar; jVar++) { - StiffMatrix_Elem[iVar][jVar] = 0.0; - } - } - - /*--- Integration formulae from "Shape functions and points of - integration of the Résumé" by Josselin Delmas (2013) ---*/ - - /*--- Tetrahedrons. Nodes of numerical integration at 1 point (order 1). ---*/ - - if (nNodes == 4) { - nGauss = 1; - Location[0][0] = 0.25; - Location[0][1] = 0.25; - Location[0][2] = 0.25; - Weight[0] = 0.166666666666666; - } - - /*--- Pyramids. Nodes numerical integration at 5 points. ---*/ - - if (nNodes == 5) { - nGauss = 5; - Location[0][0] = 0.5; - Location[0][1] = 0.0; - Location[0][2] = 0.1531754163448146; - Weight[0] = 0.133333333333333; - Location[1][0] = 0.0; - Location[1][1] = 0.5; - Location[1][2] = 0.1531754163448146; - Weight[1] = 0.133333333333333; - Location[2][0] = -0.5; - Location[2][1] = 0.0; - Location[2][2] = 0.1531754163448146; - Weight[2] = 0.133333333333333; - Location[3][0] = 0.0; - Location[3][1] = -0.5; - Location[3][2] = 0.1531754163448146; - Weight[3] = 0.133333333333333; - Location[4][0] = 0.0; - Location[4][1] = 0.0; - Location[4][2] = 0.6372983346207416; - Weight[4] = 0.133333333333333; - } - - /*--- Prism. Nodes of numerical integration at 6 points (order 3 in Xi, order 2 in Eta and Mu ). ---*/ - - if (nNodes == 6) { - nGauss = 6; - Location[0][0] = -0.577350269189626; - Location[0][1] = 0.166666666666667; - Location[0][2] = 0.166666666666667; - Weight[0] = 0.166666666666667; - Location[1][0] = -0.577350269189626; - Location[1][1] = 0.666666666666667; - Location[1][2] = 0.166666666666667; - Weight[1] = 0.166666666666667; - Location[2][0] = -0.577350269189626; - Location[2][1] = 0.166666666666667; - Location[2][2] = 0.666666666666667; - Weight[2] = 0.166666666666667; - Location[3][0] = 0.577350269189626; - Location[3][1] = 0.166666666666667; - Location[3][2] = 0.166666666666667; - Weight[3] = 0.166666666666667; - Location[4][0] = 0.577350269189626; - Location[4][1] = 0.666666666666667; - Location[4][2] = 0.166666666666667; - Weight[4] = 0.166666666666667; - Location[5][0] = 0.577350269189626; - Location[5][1] = 0.166666666666667; - Location[5][2] = 0.666666666666667; - Weight[5] = 0.166666666666667; - } - - /*--- Hexahedrons. Nodes of numerical integration at 6 points (order 3). ---*/ - - if (nNodes == 8) { - nGauss = 8; - Location[0][0] = -0.577350269189626; - Location[0][1] = -0.577350269189626; - Location[0][2] = -0.577350269189626; - Weight[0] = 1.0; - Location[1][0] = -0.577350269189626; - Location[1][1] = -0.577350269189626; - Location[1][2] = 0.577350269189626; - Weight[1] = 1.0; - Location[2][0] = -0.577350269189626; - Location[2][1] = 0.577350269189626; - Location[2][2] = -0.577350269189626; - Weight[2] = 1.0; - Location[3][0] = -0.577350269189626; - Location[3][1] = 0.577350269189626; - Location[3][2] = 0.577350269189626; - Weight[3] = 1.0; - Location[4][0] = 0.577350269189626; - Location[4][1] = -0.577350269189626; - Location[4][2] = -0.577350269189626; - Weight[4] = 1.0; - Location[5][0] = 0.577350269189626; - Location[5][1] = -0.577350269189626; - Location[5][2] = 0.577350269189626; - Weight[5] = 1.0; - Location[6][0] = 0.577350269189626; - Location[6][1] = 0.577350269189626; - Location[6][2] = -0.577350269189626; - Weight[6] = 1.0; - Location[7][0] = 0.577350269189626; - Location[7][1] = 0.577350269189626; - Location[7][2] = 0.577350269189626; - Weight[7] = 1.0; - } - - for (iGauss = 0; iGauss < nGauss; iGauss++) { - Xi = Location[iGauss][0]; - Eta = Location[iGauss][1]; - Zeta = Location[iGauss][2]; - - if (nNodes == 4) Det = ShapeFunc_Tetra(Xi, Eta, Zeta, CoordCorners, DShapeFunction); - if (nNodes == 5) Det = ShapeFunc_Pyram(Xi, Eta, Zeta, CoordCorners, DShapeFunction); - if (nNodes == 6) Det = ShapeFunc_Prism(Xi, Eta, Zeta, CoordCorners, DShapeFunction); - if (nNodes == 8) Det = ShapeFunc_Hexa(Xi, Eta, Zeta, CoordCorners, DShapeFunction); - - /*--- Compute the B Matrix ---*/ - - for (iVar = 0; iVar < 6; iVar++) - for (jVar = 0; jVar < nNodes * nVar; jVar++) B_Matrix[iVar][jVar] = 0.0; - - for (iNode = 0; iNode < nNodes; iNode++) { - B_Matrix[0][0 + iNode * nVar] = DShapeFunction[iNode][0]; - B_Matrix[1][1 + iNode * nVar] = DShapeFunction[iNode][1]; - B_Matrix[2][2 + iNode * nVar] = DShapeFunction[iNode][2]; - - B_Matrix[3][0 + iNode * nVar] = DShapeFunction[iNode][1]; - B_Matrix[3][1 + iNode * nVar] = DShapeFunction[iNode][0]; - - B_Matrix[4][1 + iNode * nVar] = DShapeFunction[iNode][2]; - B_Matrix[4][2 + iNode * nVar] = DShapeFunction[iNode][1]; - - B_Matrix[5][0 + iNode * nVar] = DShapeFunction[iNode][2]; - B_Matrix[5][2 + iNode * nVar] = DShapeFunction[iNode][0]; - } - - /*--- Impose a type of stiffness for each element ---*/ - - switch (config->GetDeform_Stiffness_Type()) { - case INVERSE_VOLUME: - E = 1.0 / ElemVolume; - break; - case SOLID_WALL_DISTANCE: - E = 1.0 / ElemDistance; - break; - case CONSTANT_STIFFNESS: - E = 1.0 / EPS; - break; - } - - Nu = config->GetDeform_Coeff(); - Mu = E / (2.0 * (1.0 + Nu)); - Lambda = Nu * E / ((1.0 + Nu) * (1.0 - 2.0 * Nu)); - - /*--- Compute the D Matrix (for plane strain and 3-D)---*/ - - D_Matrix[0][0] = Lambda + 2.0 * Mu; - D_Matrix[0][1] = Lambda; - D_Matrix[0][2] = Lambda; - D_Matrix[1][0] = Lambda; - D_Matrix[1][1] = Lambda + 2.0 * Mu; - D_Matrix[1][2] = Lambda; - D_Matrix[2][0] = Lambda; - D_Matrix[2][1] = Lambda; - D_Matrix[2][2] = Lambda + 2.0 * Mu; - D_Matrix[3][3] = Mu; - D_Matrix[4][4] = Mu; - D_Matrix[5][5] = Mu; - - /*--- Compute the BT.D Matrix ---*/ - - for (iVar = 0; iVar < nNodes * nVar; iVar++) { - for (jVar = 0; jVar < 6; jVar++) { - Aux_Matrix[iVar][jVar] = 0.0; - for (kVar = 0; kVar < 6; kVar++) Aux_Matrix[iVar][jVar] += B_Matrix[kVar][iVar] * D_Matrix[kVar][jVar]; - } - } - - /*--- Compute the BT.D.B Matrix (stiffness matrix), and add to the original - matrix using Gauss integration ---*/ - - for (iVar = 0; iVar < nNodes * nVar; iVar++) { - for (jVar = 0; jVar < nNodes * nVar; jVar++) { - for (kVar = 0; kVar < 6; kVar++) { - StiffMatrix_Elem[iVar][jVar] += Weight[iGauss] * Aux_Matrix[iVar][kVar] * B_Matrix[kVar][jVar] * fabs(Det); - } - } - } - } -} - -void CVolumetricMovement::AddFEA_StiffMatrix(CGeometry* geometry, su2double** StiffMatrix_Elem, - unsigned long PointCorners[8], unsigned short nNodes) { - unsigned short iVar, jVar, iDim, jDim; - - unsigned short nVar = geometry->GetnDim(); - - su2double** StiffMatrix_Node; - StiffMatrix_Node = new su2double*[nVar]; - for (iVar = 0; iVar < nVar; iVar++) StiffMatrix_Node[iVar] = new su2double[nVar]; - - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) StiffMatrix_Node[iVar][jVar] = 0.0; - - /*--- Transform the stiffness matrix for the hexahedral element into the - contributions for the individual nodes relative to each other. ---*/ - - for (iVar = 0; iVar < nNodes; iVar++) { - for (jVar = 0; jVar < nNodes; jVar++) { - for (iDim = 0; iDim < nVar; iDim++) { - for (jDim = 0; jDim < nVar; jDim++) { - StiffMatrix_Node[iDim][jDim] = StiffMatrix_Elem[(iVar * nVar) + iDim][(jVar * nVar) + jDim]; - } - } - - StiffMatrix.AddBlock(PointCorners[iVar], PointCorners[jVar], StiffMatrix_Node); - } - } - - /*--- Deallocate memory and exit ---*/ - - for (iVar = 0; iVar < nVar; iVar++) delete[] StiffMatrix_Node[iVar]; - delete[] StiffMatrix_Node; -} - -void CVolumetricMovement::SetBoundaryDisplacements(CGeometry* geometry, CConfig* config) { - unsigned short iDim, nDim = geometry->GetnDim(), iMarker, axis = 0; - unsigned long iPoint, total_index, iVertex; - su2double *VarCoord, MeanCoord[3] = {0.0, 0.0, 0.0}, VarIncrement = 1.0; - - /*--- Get the SU2 module. SU2_CFD will use this routine for dynamically - deforming meshes (MARKER_MOVING), while SU2_DEF will use it for deforming - meshes after imposing design variable surface deformations (DV_MARKER). ---*/ - - SU2_COMPONENT Kind_SU2 = config->GetKind_SU2(); - - /*--- If requested (no by default) impose the surface deflections in - increments and solve the grid deformation equations iteratively with - successive small deformations. ---*/ - - VarIncrement = 1.0 / ((su2double)config->GetGridDef_Nonlinear_Iter()); - - /*--- As initialization, set to zero displacements of all the surfaces except the symmetry - plane (which is treated specially, see below), internal and the send-receive boundaries ---*/ - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if (((config->GetMarker_All_KindBC(iMarker) != SYMMETRY_PLANE) && - (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) && - (config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY))) { - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - LinSysRes[total_index] = 0.0; - LinSysSol[total_index] = 0.0; - StiffMatrix.DeleteValsRowi(total_index); - } - } - } - } - - /*--- Set the known displacements, note that some points of the moving surfaces - could be on on the symmetry plane, we should specify DeleteValsRowi again (just in case) ---*/ - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if (((config->GetMarker_All_Moving(iMarker) == YES) && (Kind_SU2 == SU2_COMPONENT::SU2_CFD)) || - ((config->GetMarker_All_DV(iMarker) == YES) && (Kind_SU2 == SU2_COMPONENT::SU2_DEF)) || - ((config->GetDirectDiff() == D_DESIGN) && (Kind_SU2 == SU2_COMPONENT::SU2_CFD) && - (config->GetMarker_All_DV(iMarker) == YES)) || - ((config->GetMarker_All_DV(iMarker) == YES) && (Kind_SU2 == SU2_COMPONENT::SU2_DOT))) { - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - VarCoord = geometry->vertex[iMarker][iVertex]->GetVarCoord(); - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - LinSysRes[total_index] = SU2_TYPE::GetValue(VarCoord[iDim] * VarIncrement); - LinSysSol[total_index] = SU2_TYPE::GetValue(VarCoord[iDim] * VarIncrement); - StiffMatrix.DeleteValsRowi(total_index); - } - } - } - } - - /*--- Set to zero displacements of the normal component for the symmetry plane condition ---*/ - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE)) { - su2double* Coord_0 = nullptr; - - for (iDim = 0; iDim < nDim; iDim++) MeanCoord[iDim] = 0.0; - - /*--- Store the coord of the first point to help identify the axis. ---*/ - - iPoint = geometry->vertex[iMarker][0]->GetNode(); - Coord_0 = geometry->nodes->GetCoord(iPoint); - - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - VarCoord = geometry->nodes->GetCoord(iPoint); - for (iDim = 0; iDim < nDim; iDim++) - MeanCoord[iDim] += (VarCoord[iDim] - Coord_0[iDim]) * (VarCoord[iDim] - Coord_0[iDim]); - } - for (iDim = 0; iDim < nDim; iDim++) MeanCoord[iDim] = sqrt(MeanCoord[iDim]); - if (nDim == 3) { - if ((MeanCoord[0] <= MeanCoord[1]) && (MeanCoord[0] <= MeanCoord[2])) axis = 0; - if ((MeanCoord[1] <= MeanCoord[0]) && (MeanCoord[1] <= MeanCoord[2])) axis = 1; - if ((MeanCoord[2] <= MeanCoord[0]) && (MeanCoord[2] <= MeanCoord[1])) axis = 2; - } else { - if ((MeanCoord[0] <= MeanCoord[1])) axis = 0; - if ((MeanCoord[1] <= MeanCoord[0])) axis = 1; - } - - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - total_index = iPoint * nDim + axis; - LinSysRes[total_index] = 0.0; - LinSysSol[total_index] = 0.0; - StiffMatrix.DeleteValsRowi(total_index); - } - } - } - - /*--- Don't move the nearfield plane ---*/ - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if (config->GetMarker_All_KindBC(iMarker) == NEARFIELD_BOUNDARY) { - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - LinSysRes[total_index] = 0.0; - LinSysSol[total_index] = 0.0; - StiffMatrix.DeleteValsRowi(total_index); - } - } - } - } - - /*--- Move the FSI interfaces ---*/ - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if ((config->GetMarker_All_ZoneInterface(iMarker) == YES) && (Kind_SU2 == SU2_COMPONENT::SU2_CFD)) { - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - VarCoord = geometry->vertex[iMarker][iVertex]->GetVarCoord(); - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - LinSysRes[total_index] = SU2_TYPE::GetValue(VarCoord[iDim] * VarIncrement); - LinSysSol[total_index] = SU2_TYPE::GetValue(VarCoord[iDim] * VarIncrement); - StiffMatrix.DeleteValsRowi(total_index); - } - } - } - } -} - -void CVolumetricMovement::SetBoundaryDerivatives(CGeometry* geometry, CConfig* config, - bool ForwardProjectionDerivative) { - unsigned short iDim, iMarker; - unsigned long iPoint, total_index, iVertex; - - su2double* VarCoord; - SU2_COMPONENT Kind_SU2 = config->GetKind_SU2(); - if ((config->GetDirectDiff() == D_DESIGN) && (Kind_SU2 == SU2_COMPONENT::SU2_CFD)) { - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if ((config->GetMarker_All_DV(iMarker) == YES)) { - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - VarCoord = geometry->vertex[iMarker][iVertex]->GetVarCoord(); - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - LinSysRes[total_index] = SU2_TYPE::GetDerivative(VarCoord[iDim]); - LinSysSol[total_index] = SU2_TYPE::GetDerivative(VarCoord[iDim]); - } - } - } - } - if (LinSysRes.norm() == 0.0) cout << "Warning: Derivatives are zero!" << endl; - } else if ((Kind_SU2 == SU2_COMPONENT::SU2_DOT) && !ForwardProjectionDerivative) { - for (iPoint = 0; iPoint < nPoint; iPoint++) { - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - LinSysRes[total_index] = SU2_TYPE::GetValue(geometry->GetSensitivity(iPoint, iDim)); - LinSysSol[total_index] = SU2_TYPE::GetValue(geometry->GetSensitivity(iPoint, iDim)); - } - } - } else if (config->GetSmoothGradient() && ForwardProjectionDerivative) { - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if ((config->GetMarker_All_DV(iMarker) == YES)) { - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - LinSysRes[total_index] = SU2_TYPE::GetValue(geometry->GetSensitivity(iPoint, iDim)); - LinSysSol[total_index] = SU2_TYPE::GetValue(geometry->GetSensitivity(iPoint, iDim)); - } - } - } - } - if (LinSysRes.norm() == 0.0) cout << "Warning: Derivatives are zero!" << endl; - } -} - -void CVolumetricMovement::UpdateGridCoord_Derivatives(CGeometry* geometry, CConfig* config, - bool ForwardProjectionDerivative) { - unsigned short iDim, iMarker; - unsigned long iPoint, total_index, iVertex; - auto* new_coord = new su2double[3]; - - SU2_COMPONENT Kind_SU2 = config->GetKind_SU2(); - - /*--- Update derivatives of the grid coordinates using the solution of the linear system - after grid deformation (LinSysSol contains the derivatives of the x, y, z displacements). ---*/ - if ((config->GetDirectDiff() == D_DESIGN) && (Kind_SU2 == SU2_COMPONENT::SU2_CFD)) { - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - new_coord[0] = 0.0; - new_coord[1] = 0.0; - new_coord[2] = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - new_coord[iDim] = geometry->nodes->GetCoord(iPoint, iDim); - SU2_TYPE::SetDerivative(new_coord[iDim], SU2_TYPE::GetValue(LinSysSol[total_index])); - } - geometry->nodes->SetCoord(iPoint, new_coord); - } - } else if ((Kind_SU2 == SU2_COMPONENT::SU2_DOT) && !ForwardProjectionDerivative) { - // need to reset here, since we read out the whole vector, but are only interested in boundary derivatives. - if (config->GetSmoothGradient()) { - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - geometry->SetSensitivity(iPoint, iDim, 0.0); - } - } - } - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if (config->GetSolid_Wall(iMarker) || (config->GetMarker_All_DV(iMarker) == YES)) { - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - if (geometry->nodes->GetDomain(iPoint)) { - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - geometry->SetSensitivity(iPoint, iDim, LinSysSol[total_index]); - } - } - } - } - } - } else if (config->GetSmoothGradient() && ForwardProjectionDerivative) { - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - geometry->SetSensitivity(iPoint, iDim, LinSysSol[total_index]); - } - } - } - - delete[] new_coord; -} - -void CVolumetricMovement::SetDomainDisplacements(CGeometry* geometry, CConfig* config) { - unsigned short iDim, nDim = geometry->GetnDim(); - unsigned long iPoint, total_index; - - if (config->GetHold_GridFixed()) { - auto MinCoordValues = config->GetHold_GridFixed_Coord(); - auto MaxCoordValues = &config->GetHold_GridFixed_Coord()[3]; - - /*--- Set to zero displacements of all the points that are not going to be moved - except the surfaces ---*/ - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - auto Coord = geometry->nodes->GetCoord(iPoint); - for (iDim = 0; iDim < nDim; iDim++) { - if ((Coord[iDim] < MinCoordValues[iDim]) || (Coord[iDim] > MaxCoordValues[iDim])) { - total_index = iPoint * nDim + iDim; - LinSysRes[total_index] = 0.0; - LinSysSol[total_index] = 0.0; - StiffMatrix.DeleteValsRowi(total_index); - } - } - } - } - - /*--- Don't move the volume grid outside the limits based - on the distance to the solid surface ---*/ - - if (config->GetDeform_Limit() < 1E6) { - for (iPoint = 0; iPoint < nPoint; iPoint++) { - if (geometry->nodes->GetWall_Distance(iPoint) >= config->GetDeform_Limit()) { - for (iDim = 0; iDim < nDim; iDim++) { - total_index = iPoint * nDim + iDim; - LinSysRes[total_index] = 0.0; - LinSysSol[total_index] = 0.0; - StiffMatrix.DeleteValsRowi(total_index); - } - } - } - } -} - void CVolumetricMovement::Rigid_Rotation(CGeometry* geometry, CConfig* config, unsigned short iZone, - unsigned long iter) { + unsigned long iter) { /*--- Local variables ---*/ unsigned short iDim, nDim; unsigned long iPoint; @@ -2067,7 +630,7 @@ void CVolumetricMovement::Rigid_Rotation(CGeometry* geometry, CConfig* config, u } void CVolumetricMovement::Rigid_Pitching(CGeometry* geometry, CConfig* config, unsigned short iZone, - unsigned long iter) { + unsigned long iter) { /*--- Local variables ---*/ su2double r[3] = {0.0, 0.0, 0.0}, rotCoord[3] = {0.0, 0.0, 0.0}, *Coord, Center[3] = {0.0, 0.0, 0.0}, Omega[3] = {0.0, 0.0, 0.0}, Ampl[3] = {0.0, 0.0, 0.0}, Phase[3] = {0.0, 0.0, 0.0}; @@ -2213,7 +776,7 @@ void CVolumetricMovement::Rigid_Pitching(CGeometry* geometry, CConfig* config, u } void CVolumetricMovement::Rigid_Plunging(CGeometry* geometry, CConfig* config, unsigned short iZone, - unsigned long iter) { + unsigned long iter) { /*--- Local variables ---*/ su2double deltaX[3], newCoord[3] = {0.0, 0.0, 0.0}, Center[3], *Coord, Omega[3], Ampl[3], Lref; su2double *GridVel, newGridVel[3] = {0.0, 0.0, 0.0}, xDot[3]; @@ -2603,3 +1166,4 @@ void CVolumetricMovement::SetVolume_Rotation(CGeometry* geometry, CConfig* confi /*--- After moving all nodes, update geometry class ---*/ if (UpdateGeo) UpdateDualGrid(geometry, config); } + diff --git a/Common/src/grid_movement/CVolumetricMovementFactory.cpp b/Common/src/grid_movement/CVolumetricMovementFactory.cpp new file mode 100644 index 000000000000..62a346aa8a8b --- /dev/null +++ b/Common/src/grid_movement/CVolumetricMovementFactory.cpp @@ -0,0 +1,50 @@ +/*! + * \file CVolumetricMovementFactory.cpp + * \brief Factory to generate volumetric mover objects. + * \version 8.0.1 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/CConfig.hpp" + +#include "../../include/grid_movement/CVolumetricMovementFactory.hpp" +#include "../../include/grid_movement/CRadialBasisFunctionInterpolation.hpp" +#include "../../include/grid_movement/CLinearElasticity.hpp" + + +namespace CVolumetricMovementFactory{ + + CVolumetricMovement* CreateCVolumetricMovement(CGeometry* geometry, CConfig* config){ + const auto type = config->GetDeform_Kind(); + CVolumetricMovement* VolumetricMovement = nullptr; + switch(type){ + case DEFORM_KIND::RBF: + VolumetricMovement = new CRadialBasisFunctionInterpolation(geometry, config); + break; + case DEFORM_KIND::ELASTIC: + VolumetricMovement = new CLinearElasticity(geometry, config); + + } + return VolumetricMovement; + } + + } // namespace CVolumetricMovemementFactory diff --git a/Common/src/grid_movement/meson.build b/Common/src/grid_movement/meson.build index e543d642a944..9def44494839 100644 --- a/Common/src/grid_movement/meson.build +++ b/Common/src/grid_movement/meson.build @@ -4,4 +4,8 @@ common_src += files(['CGridMovement.cpp', 'CBezierBlending.cpp', 'CFreeFormDefBox.cpp', 'CVolumetricMovement.cpp', - 'CSurfaceMovement.cpp']) + 'CSurfaceMovement.cpp', + 'CVolumetricMovementFactory.cpp', + 'CLinearElasticity.cpp', + 'CRadialBasisFunctionInterpolation.cpp', + 'CRadialBasisFunctionNode.cpp']) diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 7160cc93e6fe..70f6f6126931 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -96,6 +96,8 @@ #include "../../../Common/include/parallelization/omp_structure.hpp" +#include "../../../Common/include/grid_movement/CVolumetricMovementFactory.hpp" + #include #ifdef VTUNEPROF @@ -2387,7 +2389,8 @@ void CDriver::PreprocessDynamicMesh(CConfig *config, CGeometry **geometry, CSolv if (!fem_solver && (config->GetGrid_Movement() || (config->GetDirectDiff() == D_DESIGN))) { if (rank == MASTER_NODE) cout << "Setting dynamic mesh structure for zone "<< iZone + 1<<"." << endl; - grid_movement = new CVolumetricMovement(geometry[MESH_0], config); + // grid_movement = new CVolumetricMovement(geometry[MESH_0], config); + grid_movement = CVolumetricMovementFactory::CreateCVolumetricMovement(geometry[MESH_0], config); surface_movement = new CSurfaceMovement(); surface_movement->CopyBoundary(geometry[MESH_0], config); diff --git a/SU2_DEF/src/drivers/CDeformationDriver.cpp b/SU2_DEF/src/drivers/CDeformationDriver.cpp index ed33ef219bc5..e620e6ac1f3f 100644 --- a/SU2_DEF/src/drivers/CDeformationDriver.cpp +++ b/SU2_DEF/src/drivers/CDeformationDriver.cpp @@ -32,6 +32,8 @@ #include "../../../SU2_CFD/include/output/CMeshOutput.hpp" #include "../../../SU2_CFD/include/solvers/CMeshSolver.hpp" +#include "../../../Common/include/grid_movement/CVolumetricMovementFactory.hpp" + using namespace std; CDeformationDriver::CDeformationDriver(char* confFile, SU2_Comm MPICommunicator) @@ -315,8 +317,8 @@ void CDeformationDriver::DeformLegacy() { /*--- Definition of the Class for grid movement. ---*/ grid_movement[iZone] = new CVolumetricMovement*[nInst_Zone](); - grid_movement[iZone][INST_0] = - new CVolumetricMovement(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); + grid_movement[iZone][INST_0] = + CVolumetricMovementFactory::CreateCVolumetricMovement(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); /*--- Save original coordinates to be reused in convexity checking procedure. ---*/ diff --git a/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp b/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp index 6e07e6b4108e..3421bf615a71 100644 --- a/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp +++ b/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp @@ -31,7 +31,7 @@ #include "../../../Common/include/geometry/CPhysicalGeometry.hpp" #include "../../../Common/include/grid_movement/CSurfaceMovement.hpp" -#include "../../../Common/include/grid_movement/CVolumetricMovement.hpp" +#include "../../../Common/include/grid_movement/CVolumetricMovementFactory.hpp" #include "../../../SU2_CFD/include/numerics/CGradSmoothing.hpp" #include "../../../SU2_CFD/include/output/CBaselineOutput.hpp" #include "../../../SU2_CFD/include/solvers/CBaselineSolver.hpp" @@ -285,8 +285,7 @@ void CDiscAdjDeformationDriver::Preprocess() { unsigned short nInst_Zone = nInst[iZone]; grid_movement[iZone] = new CVolumetricMovement*[nInst_Zone](); - grid_movement[iZone][INST_0] = - new CVolumetricMovement(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); + grid_movement[iZone][INST_0] = CVolumetricMovementFactory::CreateCVolumetricMovement(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]); /*--- Read in sensitivities from file. ---*/ diff --git a/externals/codi b/externals/codi index c6b039e5c9ed..9ca6c3828061 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit c6b039e5c9edb7675f90ffc725f9dd8e66571264 +Subproject commit 9ca6c38280610b3ea5337ca3e5b5085ee1c66b59