diff --git a/Common/include/config_structure.hpp b/Common/include/config_structure.hpp index d94242ddaf6d..e2703516568c 100644 --- a/Common/include/config_structure.hpp +++ b/Common/include/config_structure.hpp @@ -6519,60 +6519,6 @@ class CConfig { */ string GetMarker_Analyze_TagBound(unsigned short val_marker); - /*! - * \brief Set the total number of SEND_RECEIVE periodic transformations. - * \param[in] val_index - Total number of transformations. - */ - void SetnPeriodicIndex(unsigned short val_index); - - /*! - * \brief Get the total number of SEND_RECEIVE periodic transformations. - * \return Total number of transformations. - */ - unsigned short GetnPeriodicIndex(void); - - /*! - * \brief Set the rotation center for a periodic transformation. - * \param[in] val_index - Index corresponding to the periodic transformation. - * \param[in] center - Pointer to a vector containing the coordinate of the center. - */ - void SetPeriodicCenter(unsigned short val_index, su2double* center); - - /*! - * \brief Get the rotation center for a periodic transformation. - * \param[in] val_index - Index corresponding to the periodic transformation. - * \return A vector containing coordinates of the center point. - */ - su2double* GetPeriodicCenter(unsigned short val_index); - - /*! - * \brief Set the rotation angles for a periodic transformation. - * \param[in] val_index - Index corresponding to the periodic transformation. - * \param[in] rotation - Pointer to a vector containing the rotation angles. - */ - void SetPeriodicRotation(unsigned short val_index, su2double* rotation); - - /*! - * \brief Get the rotation angles for a periodic transformation. - * \param[in] val_index - Index corresponding to the periodic transformation. - * \return A vector containing the angles of rotation. - */ - su2double* GetPeriodicRotation(unsigned short val_index); - - /*! - * \brief Set the translation vector for a periodic transformation. - * \param[in] val_index - Index corresponding to the periodic transformation. - * \param[in] translate - Pointer to a vector containing the coordinate of the center. - */ - void SetPeriodicTranslate(unsigned short val_index, su2double* translate); - - /*! - * \brief Get the translation vector for a periodic transformation. - * \param[in] val_index - Index corresponding to the periodic transformation. - * \return The translation vector. - */ - su2double* GetPeriodicTranslate(unsigned short val_index); - /*! * \brief Get the total temperature at a nacelle boundary. * \param[in] val_index - Index corresponding to the inlet boundary. diff --git a/Common/include/config_structure.inl b/Common/include/config_structure.inl index 984f03494eb4..5ed8e1f2c0e2 100644 --- a/Common/include/config_structure.inl +++ b/Common/include/config_structure.inl @@ -1656,26 +1656,6 @@ inline bool CConfig::GetWrite_Conv_FSI(void) { return Write_Conv_FSI; } inline bool CConfig::GetHold_GridFixed(void) { return Hold_GridFixed; } -inline unsigned short CConfig::GetnPeriodicIndex(void) { return nPeriodic_Index; } - -inline su2double* CConfig::GetPeriodicCenter(unsigned short val_index) { return Periodic_Center[val_index]; } - -inline void CConfig::SetPeriodicCenter(unsigned short val_index, su2double* center) { - for (unsigned short i = 0; i < 3; i++) Periodic_Center[val_index][i] = center[i]; -} - -inline su2double* CConfig::GetPeriodicRotation(unsigned short val_index) { return Periodic_Rotation[val_index]; } - -inline void CConfig::SetPeriodicRotation(unsigned short val_index, su2double* rotation) { - for (unsigned short i = 0; i < 3; i++) Periodic_Rotation[val_index][i] = rotation[i]; -} - -inline su2double* CConfig::GetPeriodicTranslate(unsigned short val_index) { return Periodic_Translate[val_index]; } - -inline void CConfig::SetPeriodicTranslate(unsigned short val_index, su2double* translate) { - for (unsigned short i = 0; i < 3; i++) Periodic_Translate[val_index][i] = translate[i]; -} - inline su2double CConfig::GetCyclic_Pitch(void) { return Cyclic_Pitch; } inline su2double CConfig::GetCollective_Pitch(void) { return Collective_Pitch; } diff --git a/Common/include/geometry_structure.hpp b/Common/include/geometry_structure.hpp index dc6f32e4478d..09e1d3859617 100644 --- a/Common/include/geometry_structure.hpp +++ b/Common/include/geometry_structure.hpp @@ -3001,51 +3001,6 @@ void SetTranslationalVelocity(CConfig *config, unsigned short val_iZone, bool pr }; -/*! - * \class CPeriodicGeometry - * \brief Class for defining a periodic boundary condition. - * \author T. Economon, F. Palacios - */ -class CPeriodicGeometry : public CGeometry { - CPrimalGrid*** newBoundPer; /*!< \brief Boundary vector for new periodic elements (primal grid information). */ - unsigned long *nNewElem_BoundPer; /*!< \brief Number of new periodic elements of the boundary. */ - -public: - - /*! - * \brief Constructor of the class. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - CPeriodicGeometry(CGeometry *geometry, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CPeriodicGeometry(void); - - /*! - * \brief Set the periodic boundaries of the grid. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void SetPeriodicBoundary(CGeometry *geometry, CConfig *config); - - /*! - * \brief Set the Tecplot file. - * \param[in] config_filename - Name of the file where the Tecplot - * information is going to be stored. - */ - void SetTecPlot(char config_filename[MAX_STRING_SIZE], bool new_file); - - /*! - * \brief Write the .su2 file. - * \param[in] config - Definition of the particular problem. - * \param[in] val_mesh_out_filename - Name of the output file. - */ - void SetMeshFile(CGeometry *geometry, CConfig *config, string val_mesh_out_filename); -}; - /*! * \struct CMultiGridQueue * \brief Class for a multigrid queue system diff --git a/Common/include/matrix_structure.hpp b/Common/include/matrix_structure.hpp index b538e0970b97..c6db8c5be5a7 100644 --- a/Common/include/matrix_structure.hpp +++ b/Common/include/matrix_structure.hpp @@ -401,23 +401,6 @@ class CSysMatrix { */ void DiagonalProduct(CSysVector & vec, unsigned long row_i); - /*! - * \brief Send receive the solution using MPI. - * \param[in] x - Solution.. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - template - void SendReceive_Solution(CSysVector & x, CGeometry *geometry, CConfig *config); - - /*! - * \brief Send receive the solution using MPI and the transposed structure of the matrix. - * \param[in] x - Solution.. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void SendReceive_SolutionTransposed(CSysVector & x, CGeometry *geometry, CConfig *config); - /*! * \brief Performs the product of i-th row of a sparse matrix by a vector. * \param[in] vec - Vector to be multiplied by the row of the sparse matrix A. diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp index 9c46ea55b736..e11007cbc71e 100644 --- a/Common/src/config_structure.cpp +++ b/Common/src/config_structure.cpp @@ -1,4 +1,4 @@ -/*! +/*! * \file config_structure.cpp * \brief Main file for managing the config file * \author F. Palacios, T. Economon, B. Tracey, H. Kline @@ -4322,7 +4322,17 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ (Kind_Solver != DISC_ADJ_RANS)) { Kind_ConductivityModel_Turb = NO_CONDUCTIVITY_TURB; } - + + /*--- Check for running SU2_MSH for periodic preprocessing, and throw + an error to report that this is no longer necessary. ---*/ + + if ((Kind_SU2 == SU2_MSH) && + (Kind_Adaptation == PERIODIC)) { + SU2_MPI::Error(string("For SU2 v7.0.0 and later, preprocessing of periodic grids by SU2_MSH\n") + + string("is no longer necessary. Please use the original mesh file (prior to SU2_MSH)\n") + + string("with the same MARKER_PERIODIC definition in the configuration file.") , CURRENT_FUNCTION); + } + } void CConfig::SetMarkers(unsigned short val_software) { @@ -7862,24 +7872,6 @@ unsigned short CConfig::GetMarker_CfgFile_EngineExhaust(string val_marker) { return kMarker_All; } -void CConfig::SetnPeriodicIndex(unsigned short val_index) { - - /*--- Store total number of transformations. ---*/ - nPeriodic_Index = val_index; - - /*--- Allocate memory for centers, angles, translations. ---*/ - Periodic_Center = new su2double*[nPeriodic_Index]; - Periodic_Rotation = new su2double*[nPeriodic_Index]; - Periodic_Translate = new su2double*[nPeriodic_Index]; - - for (unsigned long i = 0; i < nPeriodic_Index; i++) { - Periodic_Center[i] = new su2double[3]; - Periodic_Rotation[i] = new su2double[3]; - Periodic_Translate[i] = new su2double[3]; - } - -} - unsigned short CConfig::GetMarker_Moving(string val_marker) { unsigned short iMarker_Moving; diff --git a/Common/src/geometry_structure.cpp b/Common/src/geometry_structure.cpp index 34c5abd1c605..62268732743c 100644 --- a/Common/src/geometry_structure.cpp +++ b/Common/src/geometry_structure.cpp @@ -9204,7 +9204,7 @@ void CPhysicalGeometry::Read_SU2_Format_Parallel(CConfig *config, string val_mes position = text_line.find ("NPERIODIC=",0); if (position != string::npos) { - unsigned short nPeriodic, iPeriodic, iIndex; + unsigned short nPeriodic; /*--- Set bool signifying that periodic transormations were found ---*/ found_transform = true; @@ -9218,36 +9218,6 @@ void CPhysicalGeometry::Read_SU2_Format_Parallel(CConfig *config, string val_mes string("is no longer necessary. Please use the original mesh file (prior to SU2_MSH)\n") + string("with the same MARKER_PERIODIC definition in the configuration file.") , CURRENT_FUNCTION); } - config->SetnPeriodicIndex(nPeriodic); - - /*--- Store center, rotation, & translation in that order for each. ---*/ - for (iPeriodic = 0; iPeriodic < nPeriodic; iPeriodic++) { - getline (mesh_file, text_line); - position = text_line.find ("PERIODIC_INDEX=",0); - if (position != string::npos) { - text_line.erase (0,15); iIndex = atoi(text_line.c_str()); - if (iIndex != iPeriodic) { - SU2_MPI::Error("PERIODIC_INDEX out of order in SU2 file!!", CURRENT_FUNCTION); - } - } - su2double* center = new su2double[3]; - su2double* rotation = new su2double[3]; - su2double* translate = new su2double[3]; - getline (mesh_file, text_line); - istringstream cent(text_line); - cent >> center[0]; cent >> center[1]; cent >> center[2]; - config->SetPeriodicCenter(iPeriodic, center); - getline (mesh_file, text_line); - istringstream rot(text_line); - rot >> rotation[0]; rot >> rotation[1]; rot >> rotation[2]; - config->SetPeriodicRotation(iPeriodic, rotation); - getline (mesh_file, text_line); - istringstream tran(text_line); - tran >> translate[0]; tran >> translate[1]; tran >> translate[2]; - config->SetPeriodicTranslate(iPeriodic, translate); - - delete [] center; delete [] rotation; delete [] translate; - } } } @@ -15061,11 +15031,9 @@ void CPhysicalGeometry::VisualizeControlVolume(CConfig *config, unsigned short a void CPhysicalGeometry::SetMeshFile (CConfig *config, string val_mesh_out_filename) { unsigned long iElem, iPoint, iElem_Bound; unsigned short iMarker, iNodes, iDim; - unsigned short iPeriodic, nPeriodic = 0; ofstream output_file; string Grid_Marker; char *cstr; - su2double *center, *angles, *transl; cstr = new char [val_mesh_out_filename.size()+1]; strcpy (cstr, val_mesh_out_filename.c_str()); @@ -15146,29 +15114,6 @@ void CPhysicalGeometry::SetMeshFile (CConfig *config, string val_mesh_out_filena } } - /*--- Get the total number of periodic transformations ---*/ - - nPeriodic = config->GetnPeriodicIndex(); - output_file << "NPERIODIC= " << nPeriodic << endl; - - /*--- From iPeriodic obtain the iMarker ---*/ - - for (iPeriodic = 0; iPeriodic < nPeriodic; iPeriodic++) { - - /*--- Retrieve the supplied periodic information. ---*/ - - center = config->GetPeriodicCenter(iPeriodic); - angles = config->GetPeriodicRotation(iPeriodic); - transl = config->GetPeriodicTranslate(iPeriodic); - - output_file << "PERIODIC_INDEX= " << iPeriodic << endl; - output_file << center[0] << "\t" << center[1] << "\t" << center[2] << endl; - output_file << angles[0] << "\t" << angles[1] << "\t" << angles[2] << endl; - output_file << transl[0] << "\t" << transl[1] << "\t" << transl[2] << endl; - - } - - output_file.close(); } @@ -15985,10 +15930,10 @@ void CPhysicalGeometry::SetGridVelocity(CConfig *config, unsigned long iter) { void CPhysicalGeometry::Set_MPI_Coord(CConfig *config) { - unsigned short iDim, iMarker, iPeriodic_Index, MarkerS, MarkerR; + unsigned short iDim, iMarker, MarkerS, MarkerR; unsigned long iVertex, iPoint, nVertexS, nVertexR, nBufferS_Vector, nBufferR_Vector; - su2double rotMatrix[3][3], *angles, theta, cosTheta, sinTheta, phi, cosPhi, sinPhi, psi, cosPsi, sinPsi, *Buffer_Receive_Coord = NULL, *Buffer_Send_Coord = NULL, *Coord = NULL, *newCoord = NULL; - su2double *translation; + su2double rotMatrix[3][3], angles[3] = {0.0,0.0,0.0}, theta, cosTheta, sinTheta, phi, cosPhi, sinPhi, psi, cosPsi, sinPsi, *Buffer_Receive_Coord = NULL, *Buffer_Send_Coord = NULL, *Coord = NULL, *newCoord = NULL; + su2double translation[3] = {0.0,0.0,0.0}; newCoord = new su2double[nDim]; #ifdef HAVE_MPI @@ -16050,12 +15995,12 @@ void CPhysicalGeometry::Set_MPI_Coord(CConfig *config) { /*--- Find point and its type of transformation ---*/ iPoint = vertex[MarkerR][iVertex]->GetNode(); - iPeriodic_Index = vertex[MarkerR][iVertex]->GetRotation_Type(); + //iPeriodic_Index = vertex[MarkerR][iVertex]->GetRotation_Type(); /*--- Retrieve the supplied periodic information. ---*/ - angles = config->GetPeriodicRotation(iPeriodic_Index); - translation = config->GetPeriodicTranslate(iPeriodic_Index); + //angles = config->GetPeriodicRotation(iPeriodic_Index); + //translation = config->GetPeriodicTranslate(iPeriodic_Index); /*--- Store angles separately for clarity. ---*/ @@ -16118,9 +16063,9 @@ void CPhysicalGeometry::Set_MPI_Coord(CConfig *config) { void CPhysicalGeometry::Set_MPI_GridVel(CConfig *config) { - unsigned short iDim, iMarker, iPeriodic_Index, MarkerS, MarkerR; + unsigned short iDim, iMarker, MarkerS, MarkerR; unsigned long iVertex, iPoint, nVertexS, nVertexR, nBufferS_Vector, nBufferR_Vector; - su2double rotMatrix[3][3], *angles, theta, cosTheta, sinTheta, phi, cosPhi, sinPhi, psi, cosPsi, sinPsi, *Buffer_Receive_GridVel = NULL, *Buffer_Send_GridVel = NULL, *GridVel = NULL, *newGridVel = NULL; + su2double rotMatrix[3][3], angles[3] = {0.0,0.0,0.0}, theta, cosTheta, sinTheta, phi, cosPhi, sinPhi, psi, cosPsi, sinPsi, *Buffer_Receive_GridVel = NULL, *Buffer_Send_GridVel = NULL, *GridVel = NULL, *newGridVel = NULL; newGridVel = new su2double[nDim]; @@ -16183,11 +16128,11 @@ void CPhysicalGeometry::Set_MPI_GridVel(CConfig *config) { /*--- Find point and its type of transformation ---*/ iPoint = vertex[MarkerR][iVertex]->GetNode(); - iPeriodic_Index = vertex[MarkerR][iVertex]->GetRotation_Type(); + //iPeriodic_Index = vertex[MarkerR][iVertex]->GetRotation_Type(); /*--- Retrieve the supplied periodic information. ---*/ - angles = config->GetPeriodicRotation(iPeriodic_Index); + //angles = config->GetPeriodicRotation(iPeriodic_Index); /*--- Store angles separately for clarity. ---*/ theta = angles[0]; phi = angles[1]; psi = angles[2]; @@ -16247,9 +16192,9 @@ void CPhysicalGeometry::Set_MPI_GridVel(CConfig *config) { void CPhysicalGeometry::Set_MPI_OldCoord(CConfig *config) { - unsigned short iDim, iMarker, iPeriodic_Index, MarkerS, MarkerR; + unsigned short iDim, iMarker, MarkerS, MarkerR; unsigned long iVertex, iPoint, nVertexS, nVertexR, nBufferS_Vector, nBufferR_Vector; - su2double rotMatrix[3][3], *angles, theta, cosTheta, sinTheta, phi, cosPhi, sinPhi, psi, cosPsi, sinPsi; + su2double rotMatrix[3][3], angles[3] = {0.0,0.0,0.0}, theta, cosTheta, sinTheta, phi, cosPhi, sinPhi, psi, cosPsi, sinPsi; su2double *Buffer_Receive_Coord_n = NULL, *Buffer_Send_Coord_n = NULL, *Coord_n = NULL, *newCoord_n = NULL; @@ -16314,11 +16259,11 @@ void CPhysicalGeometry::Set_MPI_OldCoord(CConfig *config) { /*--- Find point and its type of transformation ---*/ iPoint = vertex[MarkerR][iVertex]->GetNode(); - iPeriodic_Index = vertex[MarkerR][iVertex]->GetRotation_Type(); + //iPeriodic_Index = vertex[MarkerR][iVertex]->GetRotation_Type(); /*--- Retrieve the supplied periodic information. ---*/ - angles = config->GetPeriodicRotation(iPeriodic_Index); + //angles = config->GetPeriodicRotation(iPeriodic_Index); /*--- Store angles separately for clarity. ---*/ @@ -16439,11 +16384,11 @@ void CPhysicalGeometry::Set_MPI_OldCoord(CConfig *config) { /*--- Find point and its type of transformation ---*/ iPoint = vertex[MarkerR][iVertex]->GetNode(); - iPeriodic_Index = vertex[MarkerR][iVertex]->GetRotation_Type(); + //iPeriodic_Index = vertex[MarkerR][iVertex]->GetRotation_Type(); /*--- Retrieve the supplied periodic information. ---*/ - angles = config->GetPeriodicRotation(iPeriodic_Index); + //angles = config->GetPeriodicRotation(iPeriodic_Index); /*--- Store angles separately for clarity. ---*/ @@ -20943,629 +20888,6 @@ void CMultiGridGeometry::SetGeometryPlanes(CConfig *config) { delete[] FaceArea; } -CPeriodicGeometry::CPeriodicGeometry(CGeometry *geometry, CConfig *config) { - unsigned long nElem_new, nPoint_new, jPoint, iPoint, iElem, jElem, iVertex, - nelem_triangle = 0, nelem_quad = 0, nelem_tetra = 0, nelem_hexa = 0, nelem_prism = 0, - nelem_pyramid = 0, iIndex, newElementsBound = 0; - unsigned short iMarker, nPeriodic = 0, iPeriodic; - su2double *center, *angles, rotMatrix[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}, - translation[3], *trans, theta, phi, psi, cosTheta, sinTheta, cosPhi, sinPhi, cosPsi, sinPsi, - dx, dy, dz, rotCoord[3], *Coord_i; - unsigned short nMarker_Max = config->GetnMarker_Max(); - - /*--- We only create the mirror structure for the second boundary ---*/ - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) { - /*--- Evaluate the number of periodic boundary conditions ---*/ - nPeriodic++; - } - } - bool *CreateMirror = new bool[nPeriodic+1]; - CreateMirror[0] = false; - for (iPeriodic = 1; iPeriodic <= nPeriodic; iPeriodic++) { - if (iPeriodic <= nPeriodic/2) CreateMirror[iPeriodic] = false; - else CreateMirror[iPeriodic] = true; - } - - /*--- Write the number of dimensions of the problem ---*/ - nDim = geometry->GetnDim(); - - /*--- Copy the new boundary element information from the geometry class. - Be careful, as these are pointers to vectors/objects. ---*/ - nNewElem_BoundPer = geometry->nNewElem_Bound; - newBoundPer = geometry->newBound; - - /*--- Count the number of new boundary elements. ---*/ - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) - newElementsBound += nNewElem_BoundPer[iMarker]; - - /*--- Loop over the original grid to perform the dimensionalizaton of the new vectors ---*/ - nElem_new = 0; nPoint_new = 0; - for (iPeriodic = 1; iPeriodic <= nPeriodic; iPeriodic++) { - if (CreateMirror[iPeriodic]) { - nElem_new += geometry->PeriodicElem[iPeriodic].size(); - nPoint_new += geometry->PeriodicPoint[iPeriodic][0].size(); - } - } - - cout << "Number of new points: " << nPoint_new << "." << endl; - cout << "Number of new interior elements: " << nElem_new << "." << endl; - cout << "Number of new boundary elements added to preexisting markers: " << newElementsBound << "." << endl; - - /*--- Create a copy of the original grid ---*/ - elem = new CPrimalGrid*[geometry->GetnElem() + nElem_new]; - for (iElem = 0; iElem < geometry->GetnElem(); iElem ++) { - switch(geometry->elem[iElem]->GetVTK_Type()) { - case TRIANGLE: - elem[iElem] = new CTriangle(geometry->elem[iElem]->GetNode(0), - geometry->elem[iElem]->GetNode(1), - geometry->elem[iElem]->GetNode(2), 2); - nelem_triangle++; - break; - - case QUADRILATERAL: - elem[iElem] = new CQuadrilateral(geometry->elem[iElem]->GetNode(0), - geometry->elem[iElem]->GetNode(1), - geometry->elem[iElem]->GetNode(2), - geometry->elem[iElem]->GetNode(3), 2); - nelem_quad++; - break; - - case TETRAHEDRON: - elem[iElem] = new CTetrahedron(geometry->elem[iElem]->GetNode(0), - geometry->elem[iElem]->GetNode(1), - geometry->elem[iElem]->GetNode(2), - geometry->elem[iElem]->GetNode(3)); - nelem_tetra++; - break; - - case HEXAHEDRON: - elem[iElem] = new CHexahedron(geometry->elem[iElem]->GetNode(0), - geometry->elem[iElem]->GetNode(1), - geometry->elem[iElem]->GetNode(2), - geometry->elem[iElem]->GetNode(3), - geometry->elem[iElem]->GetNode(4), - geometry->elem[iElem]->GetNode(5), - geometry->elem[iElem]->GetNode(6), - geometry->elem[iElem]->GetNode(7)); - nelem_hexa++; - break; - - case PRISM: - elem[iElem] = new CPrism(geometry->elem[iElem]->GetNode(0), - geometry->elem[iElem]->GetNode(1), - geometry->elem[iElem]->GetNode(2), - geometry->elem[iElem]->GetNode(3), - geometry->elem[iElem]->GetNode(4), - geometry->elem[iElem]->GetNode(5)); - nelem_prism++; - break; - - case PYRAMID: - elem[iElem] = new CPyramid(geometry->elem[iElem]->GetNode(0), - geometry->elem[iElem]->GetNode(1), - geometry->elem[iElem]->GetNode(2), - geometry->elem[iElem]->GetNode(3), - geometry->elem[iElem]->GetNode(4)); - nelem_pyramid++; - break; - - } - } - - /*--- Create a list with all the points and the new index ---*/ - unsigned long *Index = new unsigned long [geometry->GetnPoint()]; - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint ++) Index[iPoint] = 0; - - for (iPeriodic = 1; iPeriodic <= nPeriodic; iPeriodic++) { - if (CreateMirror[iPeriodic]) { - for (iIndex = 0; iIndex < geometry->PeriodicPoint[iPeriodic][0].size(); iIndex++) { - iPoint = geometry->PeriodicPoint[iPeriodic][0][iIndex]; - Index[iPoint] = geometry->PeriodicPoint[iPeriodic][1][iIndex]; - } - } - } - - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) - if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - jPoint = geometry->vertex[iMarker][iVertex]->GetDonorPoint(); - Index[iPoint] = jPoint; - } - - /*--- Add the new elements due to the periodic boundary condtion ---*/ - iElem = geometry->GetnElem(); - - for (iPeriodic = 1; iPeriodic <= nPeriodic; iPeriodic++) { - if (CreateMirror[iPeriodic]) { - for (iIndex = 0; iIndex < geometry->PeriodicElem[iPeriodic].size(); iIndex++) { - jElem = geometry->PeriodicElem[iPeriodic][iIndex]; - - switch(geometry->elem[jElem]->GetVTK_Type()) { - case TRIANGLE: - elem[iElem] = new CTriangle(Index[geometry->elem[jElem]->GetNode(0)], - Index[geometry->elem[jElem]->GetNode(1)], - Index[geometry->elem[jElem]->GetNode(2)], 2); - iElem++; nelem_triangle++; - break; - - case QUADRILATERAL: - elem[iElem] = new CQuadrilateral(Index[geometry->elem[jElem]->GetNode(0)], - Index[geometry->elem[jElem]->GetNode(1)], - Index[geometry->elem[jElem]->GetNode(2)], - Index[geometry->elem[jElem]->GetNode(3)], 2); - iElem++; nelem_quad++; - break; - - case TETRAHEDRON: - elem[iElem] = new CTetrahedron(Index[geometry->elem[jElem]->GetNode(0)], - Index[geometry->elem[jElem]->GetNode(1)], - Index[geometry->elem[jElem]->GetNode(2)], - Index[geometry->elem[jElem]->GetNode(3)]); - iElem++; nelem_tetra++; - break; - - case HEXAHEDRON: - elem[iElem] = new CHexahedron(Index[geometry->elem[jElem]->GetNode(0)], - Index[geometry->elem[jElem]->GetNode(1)], - Index[geometry->elem[jElem]->GetNode(2)], - Index[geometry->elem[jElem]->GetNode(3)], - Index[geometry->elem[jElem]->GetNode(4)], - Index[geometry->elem[jElem]->GetNode(5)], - Index[geometry->elem[jElem]->GetNode(6)], - Index[geometry->elem[jElem]->GetNode(7)]); - iElem++; nelem_hexa++; - break; - - case PRISM: - elem[iElem] = new CPrism(Index[geometry->elem[jElem]->GetNode(0)], - Index[geometry->elem[jElem]->GetNode(1)], - Index[geometry->elem[jElem]->GetNode(2)], - Index[geometry->elem[jElem]->GetNode(3)], - Index[geometry->elem[jElem]->GetNode(4)], - Index[geometry->elem[jElem]->GetNode(5)]); - iElem++; nelem_prism++; - break; - - case PYRAMID: - elem[iElem] = new CPyramid(Index[geometry->elem[jElem]->GetNode(0)], - Index[geometry->elem[jElem]->GetNode(1)], - Index[geometry->elem[jElem]->GetNode(2)], - Index[geometry->elem[jElem]->GetNode(3)], - Index[geometry->elem[jElem]->GetNode(4)]); - iElem++; nelem_pyramid++; - break; - - } - } - } - } - - nElem = geometry->GetnElem() + nElem_new; - - /*--- Add the old points ---*/ - - nPointNode = geometry->GetnPoint() + nPoint_new; - node = new CPoint*[geometry->GetnPoint() + nPoint_new]; - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint ++) { - if (geometry->GetnDim() == 2) - node[iPoint] = new CPoint(geometry->node[iPoint]->GetCoord(0), - geometry->node[iPoint]->GetCoord(1), iPoint, config); - if (geometry->GetnDim() == 3) - node[iPoint] = new CPoint(geometry->node[iPoint]->GetCoord(0), - geometry->node[iPoint]->GetCoord(1), - geometry->node[iPoint]->GetCoord(2), iPoint, config); - } - - /*--- Add the new points due to the periodic boundary condtion (only in the mirror part) ---*/ - for (iPeriodic = 1; iPeriodic <= nPeriodic; iPeriodic++) { - if (CreateMirror[iPeriodic]) { - for (iIndex = 0; iIndex < geometry->PeriodicPoint[iPeriodic][0].size(); iIndex++) { - - /*--- From iPeriodic obtain the iMarker ---*/ - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) - if (iPeriodic == config->GetMarker_All_PerBound(iMarker)) break; - - /*--- Retrieve the supplied periodic information. ---*/ - center = config->GetPeriodicRotCenter(config->GetMarker_All_TagBound(iMarker)); - angles = config->GetPeriodicRotAngles(config->GetMarker_All_TagBound(iMarker)); - trans = config->GetPeriodicTranslation(config->GetMarker_All_TagBound(iMarker)); - - /*--- Store center - trans as it is constant and will be added on. - Note the subtraction, as this is the inverse translation. ---*/ - translation[0] = center[0] - trans[0]; - translation[1] = center[1] - trans[1]; - translation[2] = center[2] - trans[2]; - - /*--- Store angles separately for clarity. Compute sines/cosines. - Note the negative sign, as this is the inverse rotation. ---*/ - theta = -angles[0]; - phi = -angles[1]; - psi = -angles[2]; - - cosTheta = cos(theta); cosPhi = cos(phi); cosPsi = cos(psi); - sinTheta = sin(theta); sinPhi = sin(phi); sinPsi = sin(psi); - - /*--- Compute the rotation matrix. Note that the implicit - ordering is rotation about the x-axis, y-axis, then z-axis. ---*/ - rotMatrix[0][0] = cosPhi*cosPsi; - rotMatrix[1][0] = cosPhi*sinPsi; - rotMatrix[2][0] = -sinPhi; - - rotMatrix[0][1] = sinTheta*sinPhi*cosPsi - cosTheta*sinPsi; - rotMatrix[1][1] = sinTheta*sinPhi*sinPsi + cosTheta*cosPsi; - rotMatrix[2][1] = sinTheta*cosPhi; - - rotMatrix[0][2] = cosTheta*sinPhi*cosPsi + sinTheta*sinPsi; - rotMatrix[1][2] = cosTheta*sinPhi*sinPsi - sinTheta*cosPsi; - rotMatrix[2][2] = cosTheta*cosPhi; - - /*--- Retrieve node information for this boundary point. ---*/ - iPoint = geometry->PeriodicPoint[iPeriodic][0][iIndex]; - jPoint = geometry->PeriodicPoint[iPeriodic][1][iIndex]; - Coord_i = geometry->node[iPoint]->GetCoord(); - - /*--- Get the position vector from rot center to point. ---*/ - dx = Coord_i[0] - center[0]; - dy = Coord_i[1] - center[1]; - if (nDim == 3) { - dz = Coord_i[2] - center[2]; - } else { - dz = 0.0; - } - - /*--- Compute transformed point coordinates. ---*/ - rotCoord[0] = rotMatrix[0][0]*dx + rotMatrix[0][1]*dy + rotMatrix[0][2]*dz + translation[0]; - rotCoord[1] = rotMatrix[1][0]*dx + rotMatrix[1][1]*dy + rotMatrix[1][2]*dz + translation[1]; - rotCoord[2] = rotMatrix[2][0]*dx + rotMatrix[2][1]*dy + rotMatrix[2][2]*dz + translation[2]; - - /*--- Save the new points with the new coordinates. ---*/ - if (geometry->GetnDim() == 2) - node[jPoint] = new CPoint(rotCoord[0], rotCoord[1], jPoint, config); - if (geometry->GetnDim() == 3) - node[jPoint] = new CPoint(rotCoord[0], rotCoord[1], rotCoord[2], jPoint, config); - - } - } - } - - nPoint = geometry->GetnPoint() + nPoint_new; - - /*--- Add the old boundary, reserving space for two new bc (send/recive periodic bc) ---*/ - nMarker = geometry->GetnMarker() + 2; - nElem_Bound = new unsigned long [nMarker]; - bound = new CPrimalGrid**[nMarker]; - Tag_to_Marker = new string [nMarker_Max]; - config->SetnMarker_All(nMarker); - - /*--- Copy the olf boundary ---*/ - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - - bound[iMarker] = new CPrimalGrid* [geometry->GetnElem_Bound(iMarker)]; - - for (iVertex = 0; iVertex < geometry->GetnElem_Bound(iMarker); iVertex++) { - if (geometry->bound[iMarker][iVertex]->GetVTK_Type() == LINE) - bound[iMarker][iVertex] = new CLine(geometry->bound[iMarker][iVertex]->GetNode(0), - geometry->bound[iMarker][iVertex]->GetNode(1), 2); - if (geometry->bound[iMarker][iVertex]->GetVTK_Type() == TRIANGLE) - bound[iMarker][iVertex] = new CTriangle(geometry->bound[iMarker][iVertex]->GetNode(0), - geometry->bound[iMarker][iVertex]->GetNode(1), - geometry->bound[iMarker][iVertex]->GetNode(2), 3); - if (geometry->bound[iMarker][iVertex]->GetVTK_Type() == QUADRILATERAL) - bound[iMarker][iVertex] = new CQuadrilateral(geometry->bound[iMarker][iVertex]->GetNode(0), - geometry->bound[iMarker][iVertex]->GetNode(1), - geometry->bound[iMarker][iVertex]->GetNode(2), - geometry->bound[iMarker][iVertex]->GetNode(3), 3); - } - - nElem_Bound[iMarker] = geometry->GetnElem_Bound(iMarker); - Tag_to_Marker[iMarker] = geometry->GetMarker_Tag(iMarker); - - } - - delete [] Index; - delete [] CreateMirror; - -} - -CPeriodicGeometry::~CPeriodicGeometry(void) { - unsigned long iElem_Bound; - unsigned short iMarker; - - for (iMarker = 0; iMarker < nMarker; iMarker++) { - for (iElem_Bound = 0; iElem_Bound < nElem_Bound[iMarker]; iElem_Bound++) { - if (newBoundPer[iMarker][iElem_Bound] != NULL) delete [] newBoundPer[iMarker][iElem_Bound]; - } - } - if (newBoundPer != NULL) delete[] newBoundPer; - - if (nNewElem_BoundPer != NULL) delete[] nNewElem_BoundPer; - -} - -void CPeriodicGeometry::SetPeriodicBoundary(CGeometry *geometry, CConfig *config) { - unsigned short iMarker, iPeriodic, nPeriodic = 0, iMarkerSend, iMarkerReceive; - unsigned long iVertex, Counter_Send = 0, Counter_Receive = 0, iIndex; - - /*--- Compute the number of periodic bc on the geometry ---*/ - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) - if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) - nPeriodic++; - - /*--- First compute the Send/Receive boundaries, count the number of points ---*/ - Counter_Send = 0; Counter_Receive = 0; - for (iPeriodic = 1; iPeriodic <= nPeriodic; iPeriodic++) { - if (geometry->PeriodicPoint[iPeriodic][0].size() != 0) - Counter_Send += geometry->PeriodicPoint[iPeriodic][0].size(); - if (geometry->PeriodicPoint[iPeriodic][1].size() != 0) - Counter_Receive += geometry->PeriodicPoint[iPeriodic][1].size(); - } - - /*--- Adimensionalization of the new boundaries ---*/ - iMarkerSend = nMarker - 2; iMarkerReceive = nMarker - 1; - config->SetMarker_All_SendRecv(iMarkerSend,1); - config->SetMarker_All_SendRecv(iMarkerReceive,-1); - nElem_Bound[iMarkerSend] = Counter_Send; - nElem_Bound[iMarkerReceive] = Counter_Receive; - bound[iMarkerSend] = new CPrimalGrid* [Counter_Send]; - bound[iMarkerReceive] = new CPrimalGrid* [Counter_Receive]; - - /*--- First we do the send ---*/ - iVertex = 0; - for (iPeriodic = 1; iPeriodic <= nPeriodic; iPeriodic++) - if (geometry->PeriodicPoint[iPeriodic][0].size() != 0) - for (iIndex = 0; iIndex < geometry->PeriodicPoint[iPeriodic][0].size(); iIndex++) { - bound[iMarkerSend][iVertex] = new CVertexMPI(geometry->PeriodicPoint[iPeriodic][0][iIndex], nDim); - bound[iMarkerSend][iVertex]->SetRotation_Type(iPeriodic); - iVertex++; - } - - /*--- Second we do the receive ---*/ - iVertex = 0; - for (iPeriodic = 1; iPeriodic <= nPeriodic; iPeriodic++) - if (geometry->PeriodicPoint[iPeriodic][1].size() != 0) - for (iIndex = 0; iIndex < geometry->PeriodicPoint[iPeriodic][1].size(); iIndex++) { - bound[iMarkerReceive][iVertex] = new CVertexMPI(geometry->PeriodicPoint[iPeriodic][1][iIndex], nDim); - bound[iMarkerReceive][iVertex]->SetRotation_Type(iPeriodic); - iVertex++; - } - -} - -void CPeriodicGeometry::SetMeshFile(CGeometry *geometry, CConfig *config, string val_mesh_out_filename) { - unsigned long iElem, iPoint, iElem_Bound, GhostPoints; - unsigned short iMarker, iNodes, iDim; - unsigned short iMarkerReceive, iPeriodic, nPeriodic = 0; - ofstream output_file; - string Grid_Marker; - char *cstr; - su2double *center, *angles, *transl; - - cstr = new char [val_mesh_out_filename.size()+1]; - strcpy (cstr, val_mesh_out_filename.c_str()); - - /*--- Open .su2 grid file ---*/ - output_file.precision(15); - output_file.open(cstr, ios::out); - - /*--- Ghost points, look at the nodes in the send receive ---*/ - iMarkerReceive = nMarker - 1; - GhostPoints = nElem_Bound[iMarkerReceive]; - - /*--- Change the numbering to guarantee that the all the receive - points are at the end of the file. ---*/ - std::vector receive_nodes; - std::vector send_nodes; - for (iMarker = 0; iMarker < nMarker; iMarker++) { - if (bound[iMarker][0]->GetVTK_Type() == VERTEX) { - if (config->GetMarker_All_SendRecv(iMarker) < 0) { - for (iElem_Bound = 0; iElem_Bound < nElem_Bound[iMarker]; iElem_Bound++) { - if (bound[iMarker][iElem_Bound]->GetRotation_Type() == 1) { - receive_nodes.push_back(bound[iMarker][iElem_Bound]->GetNode(0)); - } else { - send_nodes.push_back(bound[iMarker][iElem_Bound]->GetNode(0)); - } - } - } - } - } - - /*--- Build the sorted lists of node numbers with receive/send at the end - * NewSort[i] = j maps the new number (i) to the old number (j) - * ReverseSort[j] = i maps the old number (j) to the new number (i) ---*/ - std::vector NewSort; - std::vector ReverseSort; - for (iPoint = 0; iPoint < nPoint; iPoint++) { - bool isReceive = (find(receive_nodes.begin(), receive_nodes.end(), iPoint) - != receive_nodes.end()); - bool isSend = (find(send_nodes.begin(), send_nodes.end(), iPoint) - != send_nodes.end()); - if (!isSend && !isReceive) { - NewSort.push_back(iPoint); - } - } - NewSort.insert(NewSort.end(), receive_nodes.begin(), receive_nodes.end()); - NewSort.insert(NewSort.end(), send_nodes.begin(), send_nodes.end()); - - ReverseSort.resize(NewSort.size()); - for (iPoint = 0; iPoint < nPoint; iPoint++) { - unsigned long jPoint; - for (jPoint = 0; jPoint < nPoint; jPoint++) { - if (NewSort[iPoint] == jPoint) { - ReverseSort[jPoint] = iPoint; - break; - } - } - if (jPoint == nPoint) { // Loop fell through without break - SU2_MPI::Error("Remapping of periodic nodes failed.", CURRENT_FUNCTION); - } - } - - /*--- Write dimension, number of elements and number of points ---*/ - output_file << "NDIME= " << nDim << endl; - output_file << "NELEM= " << nElem << endl; - for (iElem = 0; iElem < nElem; iElem++) { - output_file << elem[iElem]->GetVTK_Type(); - for (iNodes = 0; iNodes < elem[iElem]->GetnNodes(); iNodes++) - output_file << "\t" << ReverseSort[elem[iElem]->GetNode(iNodes)]; - output_file << "\t"<GetCoord(iDim) ; - } - output_file << "\t" << iPoint << endl; - } - - output_file << "NMARK= " << nMarker << endl; - for (iMarker = 0; iMarker < nMarker; iMarker++) { - if (bound[iMarker][0]->GetVTK_Type() != VERTEX) { - - Grid_Marker = config->GetMarker_All_TagBound(iMarker); - output_file << "MARKER_TAG= " << Grid_Marker << endl; - output_file << "MARKER_ELEMS= " << nElem_Bound[iMarker] + nNewElem_BoundPer[iMarker] << endl; - - for (iElem_Bound = 0; iElem_Bound < nElem_Bound[iMarker]; iElem_Bound++) { - output_file << bound[iMarker][iElem_Bound]->GetVTK_Type() << "\t" ; - for (iNodes = 0; iNodes < bound[iMarker][iElem_Bound]->GetnNodes()-1; iNodes++) - output_file << ReverseSort[bound[iMarker][iElem_Bound]->GetNode(iNodes)] << "\t" ; - iNodes = bound[iMarker][iElem_Bound]->GetnNodes()-1; - output_file << ReverseSort[bound[iMarker][iElem_Bound]->GetNode(iNodes)] << endl; - } - - /*--- Write any new elements at the end of the list. ---*/ - if (nNewElem_BoundPer[iMarker] > 0) { - for (iElem_Bound = 0; iElem_Bound < nNewElem_BoundPer[iMarker]; iElem_Bound++) { - output_file << newBoundPer[iMarker][iElem_Bound]->GetVTK_Type() << "\t" ; - for (iNodes = 0; iNodes < newBoundPer[iMarker][iElem_Bound]->GetnNodes()-1; iNodes++) - output_file << ReverseSort[newBoundPer[iMarker][iElem_Bound]->GetNode(iNodes)] << "\t" ; - iNodes = newBoundPer[iMarker][iElem_Bound]->GetnNodes()-1; - output_file << ReverseSort[newBoundPer[iMarker][iElem_Bound]->GetNode(iNodes)] << endl; - } - } - - } - - if (bound[iMarker][0]->GetVTK_Type() == VERTEX) { - output_file << "MARKER_TAG= SEND_RECEIVE" << endl; - output_file << "MARKER_ELEMS= " << nElem_Bound[iMarker]<< endl; - if (config->GetMarker_All_SendRecv(iMarker) > 0) output_file << "SEND_TO= " << config->GetMarker_All_SendRecv(iMarker) << endl; - if (config->GetMarker_All_SendRecv(iMarker) < 0) output_file << "SEND_TO= " << config->GetMarker_All_SendRecv(iMarker) << endl; - - for (iElem_Bound = 0; iElem_Bound < nElem_Bound[iMarker]; iElem_Bound++) { - output_file << bound[iMarker][iElem_Bound]->GetVTK_Type() << "\t" << - ReverseSort[bound[iMarker][iElem_Bound]->GetNode(0)] << "\t" << - bound[iMarker][iElem_Bound]->GetRotation_Type() << endl; - } - } - } - - /*--- Compute the number of periodic bc on the geometry ---*/ - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) - if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) - nPeriodic++; - - output_file << "NPERIODIC= " << nPeriodic + 1 << endl; - - /*--- Periodic 0 correspond with no movement of the surface ---*/ - output_file << "PERIODIC_INDEX= 0" << endl; - output_file << "0.000000000000000e+00" << "\t" << "0.000000000000000e+00" << "\t" << "0.000000000000000e+00" << endl; - output_file << "0.000000000000000e+00" << "\t" << "0.000000000000000e+00" << "\t" << "0.000000000000000e+00" << endl; - output_file << "0.000000000000000e+00" << "\t" << "0.000000000000000e+00" << "\t" << "0.000000000000000e+00" << endl; - - /*--- From iPeriodic obtain the iMarker ---*/ - for (iPeriodic = 1; iPeriodic <= nPeriodic; iPeriodic++) { - for (iMarker = 0; iMarker < nMarker; iMarker++) - if (iPeriodic == config->GetMarker_All_PerBound(iMarker)) break; - - /*--- Retrieve the supplied periodic information. ---*/ - center = config->GetPeriodicRotCenter(config->GetMarker_All_TagBound(iMarker)); - angles = config->GetPeriodicRotAngles(config->GetMarker_All_TagBound(iMarker)); - transl = config->GetPeriodicTranslation(config->GetMarker_All_TagBound(iMarker)); - - output_file << "PERIODIC_INDEX= " << iPeriodic << endl; - output_file << center[0] << "\t" << center[1] << "\t" << center[2] << endl; - output_file << angles[0] << "\t" << angles[1] << "\t" << angles[2] << endl; - output_file << transl[0] << "\t" << transl[1] << "\t" << transl[2] << endl; - - } - - - output_file.close(); -} - -void CPeriodicGeometry::SetTecPlot(char mesh_filename[MAX_STRING_SIZE], bool new_file) { - - unsigned long iElem, iPoint; - unsigned short iDim; - ofstream Tecplot_File; - - Tecplot_File.open(mesh_filename, ios::out); - Tecplot_File << "TITLE= \"Visualization of the volumetric grid\"" << endl; - - if (nDim == 2) { - Tecplot_File << "VARIABLES = \"x\",\"y\" " << endl; - Tecplot_File << "ZONE NODES= "<< nPoint <<", ELEMENTS= "<< nElem <<", DATAPACKING=POINT, ZONETYPE=FEQUADRILATERAL"<< endl; - } - if (nDim == 3) { - Tecplot_File << "VARIABLES = \"x\",\"y\",\"z\" " << endl; - Tecplot_File << "ZONE NODES= "<< nPoint <<", ELEMENTS= "<< nElem <<", DATAPACKING=POINT, ZONETYPE=FEBRICK"<< endl; - } - - for (iPoint = 0; iPoint < nPoint; iPoint++) { - for (iDim = 0; iDim < nDim; iDim++) - Tecplot_File << scientific << node[iPoint]->GetCoord(iDim) << "\t"; - Tecplot_File << "\n"; - } - - for (iElem = 0; iElem < nElem; iElem++) { - if (elem[iElem]->GetVTK_Type() == TRIANGLE) { - Tecplot_File << - elem[iElem]->GetNode(0)+1 <<" "<< elem[iElem]->GetNode(1)+1 <<" "<< - elem[iElem]->GetNode(2)+1 <<" "<< elem[iElem]->GetNode(2)+1 << endl; - } - if (elem[iElem]->GetVTK_Type() == QUADRILATERAL) { - Tecplot_File << - elem[iElem]->GetNode(0)+1 <<" "<< elem[iElem]->GetNode(1)+1 <<" "<< - elem[iElem]->GetNode(2)+1 <<" "<< elem[iElem]->GetNode(3)+1 << endl; - } - if (elem[iElem]->GetVTK_Type() == TETRAHEDRON) { - Tecplot_File << - elem[iElem]->GetNode(0)+1 <<" "<< elem[iElem]->GetNode(1)+1 <<" "<< - elem[iElem]->GetNode(2)+1 <<" "<< elem[iElem]->GetNode(2)+1 <<" "<< - elem[iElem]->GetNode(3)+1 <<" "<< elem[iElem]->GetNode(3)+1 <<" "<< - elem[iElem]->GetNode(3)+1 <<" "<< elem[iElem]->GetNode(3)+1 << endl; - } - if (elem[iElem]->GetVTK_Type() == HEXAHEDRON) { - Tecplot_File << - elem[iElem]->GetNode(0)+1 <<" "<< elem[iElem]->GetNode(1)+1 <<" "<< - elem[iElem]->GetNode(2)+1 <<" "<< elem[iElem]->GetNode(3)+1 <<" "<< - elem[iElem]->GetNode(4)+1 <<" "<< elem[iElem]->GetNode(5)+1 <<" "<< - elem[iElem]->GetNode(6)+1 <<" "<< elem[iElem]->GetNode(7)+1 << endl; - } - if (elem[iElem]->GetVTK_Type() == PYRAMID) { - Tecplot_File << - elem[iElem]->GetNode(0)+1 <<" "<< elem[iElem]->GetNode(1)+1 <<" "<< - elem[iElem]->GetNode(2)+1 <<" "<< elem[iElem]->GetNode(3)+1 <<" "<< - elem[iElem]->GetNode(4)+1 <<" "<< elem[iElem]->GetNode(4)+1 <<" "<< - elem[iElem]->GetNode(4)+1 <<" "<< elem[iElem]->GetNode(4)+1 << endl; - } - if (elem[iElem]->GetVTK_Type() == PRISM) { - Tecplot_File << - elem[iElem]->GetNode(0)+1 <<" "<< elem[iElem]->GetNode(1)+1 <<" "<< - elem[iElem]->GetNode(1)+1 <<" "<< elem[iElem]->GetNode(2)+1 <<" "<< - elem[iElem]->GetNode(3)+1 <<" "<< elem[iElem]->GetNode(4)+1 <<" "<< - elem[iElem]->GetNode(4)+1 <<" "<< elem[iElem]->GetNode(5)+1 << endl; - } - } - - Tecplot_File.close(); -} - CMultiGridQueue::CMultiGridQueue(unsigned long val_npoint) { unsigned long iPoint; diff --git a/Common/src/matrix_structure.cpp b/Common/src/matrix_structure.cpp index 840da34e44ef..4c0e4312ffba 100644 --- a/Common/src/matrix_structure.cpp +++ b/Common/src/matrix_structure.cpp @@ -1169,185 +1169,6 @@ void CSysMatrix::DiagonalProduct(CSysVector & vec, unsig } -template -template -void CSysMatrix::SendReceive_Solution(CSysVector & x, CGeometry *geometry, CConfig *config) { - - unsigned short iVar, iMarker, MarkerS, MarkerR; - unsigned long iVertex, iPoint, nVertexS, nVertexR, nBufferS_Vector, nBufferR_Vector; - OtherType *Buffer_Receive = NULL, *Buffer_Send = NULL; - -#ifdef HAVE_MPI - int send_to, receive_from; - SU2_MPI::Status status; -#endif - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - - if ((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && - (config->GetMarker_All_SendRecv(iMarker) > 0)) { - - MarkerS = iMarker; MarkerR = iMarker+1; - -#ifdef HAVE_MPI - - send_to = config->GetMarker_All_SendRecv(MarkerS)-1; - receive_from = abs(config->GetMarker_All_SendRecv(MarkerR))-1; - -#endif - - nVertexS = geometry->nVertex[MarkerS]; nVertexR = geometry->nVertex[MarkerR]; - nBufferS_Vector = nVertexS*nVar; nBufferR_Vector = nVertexR*nVar; - - /*--- Allocate Receive and send buffers ---*/ - - Buffer_Receive = new OtherType [nBufferR_Vector]; - Buffer_Send = new OtherType[nBufferS_Vector]; - - /*--- Copy the solution that should be sended ---*/ - - for (iVertex = 0; iVertex < nVertexS; iVertex++) { - iPoint = geometry->vertex[MarkerS][iVertex]->GetNode(); - for (iVar = 0; iVar < nVar; iVar++) - Buffer_Send[iVertex*nVar+iVar] = x[iPoint*nVar+iVar]; - } - -#ifdef HAVE_MPI - - /*--- Send/Receive information using Sendrecv ---*/ - - SelectMPIWrapper::W::Sendrecv(Buffer_Send, nBufferS_Vector, MPI_DOUBLE, send_to, 0, - Buffer_Receive, nBufferR_Vector, MPI_DOUBLE, receive_from, 0, MPI_COMM_WORLD, &status); - -#else - - /*--- Receive information without MPI ---*/ - - for (iVertex = 0; iVertex < nVertexR; iVertex++) { - for (iVar = 0; iVar < nVar; iVar++) - Buffer_Receive[iVar*nVertexR+iVertex] = Buffer_Send[iVar*nVertexR+iVertex]; - } - -#endif - - /*--- Deallocate send buffer ---*/ - - delete [] Buffer_Send; - - /*--- Do the coordinate transformation ---*/ - - for (iVertex = 0; iVertex < nVertexR; iVertex++) { - - /*--- Find point and its type of transformation ---*/ - - iPoint = geometry->vertex[MarkerR][iVertex]->GetNode(); - - /*--- Copy transformed conserved variables back into buffer. ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - x[iPoint*nVar+iVar] = Buffer_Receive[iVertex*nVar+iVar]; - - } - - /*--- Deallocate receive buffer ---*/ - - delete [] Buffer_Receive; - - } - - } - -} - -template -void CSysMatrix::SendReceive_SolutionTransposed(CSysVector & x, CGeometry *geometry, CConfig *config) { - - unsigned short iVar, iMarker, MarkerS, MarkerR; - unsigned long iVertex, iPoint, nVertexS, nVertexR, nBufferS_Vector, nBufferR_Vector; - ScalarType *Buffer_Receive = NULL, *Buffer_Send = NULL; - -#ifdef HAVE_MPI - int send_to, receive_from; - SU2_MPI::Status status; -#endif - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - - if ((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && - (config->GetMarker_All_SendRecv(iMarker) > 0)) { - - MarkerS = iMarker + 1; MarkerR = iMarker; - -#ifdef HAVE_MPI - - receive_from = config->GetMarker_All_SendRecv(MarkerR)-1; - send_to = abs(config->GetMarker_All_SendRecv(MarkerS))-1; - -#endif - - nVertexS = geometry->nVertex[MarkerS]; nVertexR = geometry->nVertex[MarkerR]; - nBufferS_Vector = nVertexS*nVar; nBufferR_Vector = nVertexR*nVar; - - /*--- Allocate Receive and send buffers ---*/ - - Buffer_Receive = new ScalarType [nBufferR_Vector]; - Buffer_Send = new ScalarType[nBufferS_Vector]; - - /*--- Copy the solution that should be sended ---*/ - - for (iVertex = 0; iVertex < nVertexS; iVertex++) { - iPoint = geometry->vertex[MarkerS][iVertex]->GetNode(); - for (iVar = 0; iVar < nVar; iVar++) - Buffer_Send[iVertex*nVar+iVar] = x[iPoint*nVar+iVar]; - } - -#ifdef HAVE_MPI - - /*--- Send/Receive information using Sendrecv ---*/ - - SelectMPIWrapper::W::Sendrecv(Buffer_Send, nBufferS_Vector, MPI_DOUBLE, send_to, 0, - Buffer_Receive, nBufferR_Vector, MPI_DOUBLE, receive_from, 0, MPI_COMM_WORLD, &status); - -#else - - /*--- Receive information without MPI ---*/ - - for (iVertex = 0; iVertex < nVertexR; iVertex++) { - for (iVar = 0; iVar < nVar; iVar++) - Buffer_Receive[iVar*nVertexR+iVertex] = Buffer_Send[iVar*nVertexR+iVertex]; - } - -#endif - - /*--- Deallocate send buffer ---*/ - - delete [] Buffer_Send; - - /*--- Do the coordinate transformation ---*/ - - for (iVertex = 0; iVertex < nVertexR; iVertex++) { - - /*--- Find point and its type of transformation ---*/ - - iPoint = geometry->vertex[MarkerR][iVertex]->GetNode(); - - /*--- Copy transformed conserved variables back into buffer. ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - x[iPoint*nVar+iVar] += Buffer_Receive[iVertex*nVar+iVar]; - - } - - /*--- Deallocate receive buffer ---*/ - - delete [] Buffer_Receive; - - } - - } - -} - template void CSysMatrix::RowProduct(const CSysVector & vec, unsigned long row_i) { @@ -2570,7 +2391,6 @@ void CSysMatrix::ComputeResidual(const CSysVector & sol, /*--- Explicit instantiations ---*/ template class CSysMatrix; -template void CSysMatrix::SendReceive_Solution(CSysVector&, CGeometry*, CConfig*); template void CSysMatrix::InitiateComms(CSysVector&, CGeometry*, CConfig*, unsigned short); template void CSysMatrix::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short); template class CSysMatrixVectorProduct; @@ -2582,8 +2402,6 @@ template class CLineletPreconditioner; #ifdef CODI_REVERSE_TYPE template class CSysMatrix; -template void CSysMatrix::SendReceive_Solution(CSysVector&, CGeometry*, CConfig*); -template void CSysMatrix::SendReceive_Solution(CSysVector&, CGeometry*, CConfig*); template void CSysMatrix::InitiateComms(CSysVector&, CGeometry*, CConfig*, unsigned short); template void CSysMatrix::InitiateComms(CSysVector&, CGeometry*, CConfig*, unsigned short); template void CSysMatrix::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short); diff --git a/SU2_CFD/src/output_su2.cpp b/SU2_CFD/src/output_su2.cpp index f3994e9d1900..fc30602387f4 100644 --- a/SU2_CFD/src/output_su2.cpp +++ b/SU2_CFD/src/output_su2.cpp @@ -40,9 +40,8 @@ void COutput::SetSU2_MeshASCII(CConfig *config, CGeometry *geometry, unsigned short val_iZone, ofstream& output_file) { unsigned long iElem, iPoint, iElem_Bound, nElem_Bound_, vnodes_edge[2], vnodes_triangle[3], vnodes_quad[4], iNode, nElem; - unsigned short iMarker, iDim, nDim = geometry->GetnDim(), iChar, iPeriodic, nPeriodic = 0, VTK_Type, nMarker_; + unsigned short iMarker, iDim, nDim = geometry->GetnDim(), iChar, VTK_Type, nMarker_; short SendTo; - su2double *center, *angles, *transl; ifstream input_file; string Grid_Marker, text_line, Marker_Tag, str; string::size_type position; @@ -221,28 +220,6 @@ void COutput::SetSU2_MeshASCII(CConfig *config, CGeometry *geometry, unsigned sh remove(str.c_str()); - /*--- Get the total number of periodic transformations ---*/ - - nPeriodic = config->GetnPeriodicIndex(); - output_file << "NPERIODIC= " << nPeriodic << endl; - - /*--- From iPeriodic obtain the iMarker ---*/ - - for (iPeriodic = 0; iPeriodic < nPeriodic; iPeriodic++) { - - /*--- Retrieve the supplied periodic information. ---*/ - - center = config->GetPeriodicCenter(iPeriodic); - angles = config->GetPeriodicRotation(iPeriodic); - transl = config->GetPeriodicTranslate(iPeriodic); - - output_file << "PERIODIC_INDEX= " << iPeriodic << endl; - output_file << center[0] << "\t" << center[1] << "\t" << center[2] << endl; - output_file << angles[0] << "\t" << angles[1] << "\t" << angles[2] << endl; - output_file << transl[0] << "\t" << transl[1] << "\t" << transl[2] << endl; - - } - } void COutput::SetSU2_MeshBinary(CConfig *config, CGeometry *geometry) { } diff --git a/SU2_CFD/src/solver_structure.cpp b/SU2_CFD/src/solver_structure.cpp index 3213ac1c8c1b..8a4469597e9c 100644 --- a/SU2_CFD/src/solver_structure.cpp +++ b/SU2_CFD/src/solver_structure.cpp @@ -3550,91 +3550,6 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, CConfig *config) { } -void CSolver::SetPressureLaplacian(CGeometry *geometry, CConfig *config, su2double *PressureLaplacian) { - - unsigned long Point = 0, iPoint = 0, jPoint = 0, iEdge, iVertex; - unsigned short iMarker, iVar; - su2double DualArea, Partial_Res, *Normal; - su2double **UxVar_Gradient, **UyVar_Gradient; - - UxVar_Gradient = new su2double* [geometry->GetnPoint()]; - UyVar_Gradient = new su2double* [geometry->GetnPoint()]; - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - UxVar_Gradient[iPoint] = new su2double [2]; - UyVar_Gradient[iPoint] = new su2double [2]; - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) - for (iVar = 0; iVar < 2; iVar++) { - UxVar_Gradient[iPoint][iVar] = 0.0; - UyVar_Gradient[iPoint][iVar] = 0.0; - } - - /*--- Loop interior edges ---*/ - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Normal = geometry->edge[iEdge]->GetNormal(); - - Partial_Res = 0.5 * ( node[iPoint]->GetSolution(1) + node[jPoint]->GetSolution(1)) * Normal[0]; - UxVar_Gradient[iPoint][0] += Partial_Res; - UxVar_Gradient[jPoint][0] -= Partial_Res; - - Partial_Res = 0.5 * ( node[iPoint]->GetSolution(1) + node[jPoint]->GetSolution(1)) * Normal[1]; - UxVar_Gradient[iPoint][1] += Partial_Res; - UxVar_Gradient[jPoint][1] -= Partial_Res; - - Partial_Res = 0.5 * ( node[iPoint]->GetSolution(2) + node[jPoint]->GetSolution(2)) * Normal[0]; - UyVar_Gradient[iPoint][0] += Partial_Res; - UyVar_Gradient[jPoint][0] -= Partial_Res; - - Partial_Res = 0.5 * ( node[iPoint]->GetSolution(2) + node[jPoint]->GetSolution(2)) * Normal[1]; - UyVar_Gradient[iPoint][1] += Partial_Res; - UyVar_Gradient[jPoint][1] -= Partial_Res; - - } - - /*--- Loop boundary edges ---*/ - - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) - if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && - (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - Point = geometry->vertex[iMarker][iVertex]->GetNode(); - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - - Partial_Res = node[Point]->GetSolution(1) * Normal[0]; - UxVar_Gradient[Point][0] -= Partial_Res; - - Partial_Res = node[Point]->GetSolution(1) * Normal[1]; - UxVar_Gradient[Point][1] -= Partial_Res; - - Partial_Res = node[Point]->GetSolution(2) * Normal[0]; - UyVar_Gradient[Point][0] -= Partial_Res; - - Partial_Res = node[Point]->GetSolution(2) * Normal[1]; - UyVar_Gradient[Point][1] -= Partial_Res; - } - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - DualArea = geometry->node[iPoint]->GetVolume(); - PressureLaplacian[iPoint] = (UxVar_Gradient[iPoint][0]*UxVar_Gradient[iPoint][0] + UyVar_Gradient[iPoint][1]*UyVar_Gradient[iPoint][1] + - UxVar_Gradient[iPoint][1]*UyVar_Gradient[iPoint][0] + UxVar_Gradient[iPoint][0]*UyVar_Gradient[iPoint][1])/DualArea ; - - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - delete[] UxVar_Gradient[iPoint]; - delete[] UyVar_Gradient[iPoint]; - } - - delete[] UxVar_Gradient; - delete[] UyVar_Gradient; - -} - void CSolver::Gauss_Elimination(su2double** A, su2double* rhs, unsigned short nVar) { short iVar, jVar, kVar; diff --git a/SU2_MSH/src/SU2_MSH.cpp b/SU2_MSH/src/SU2_MSH.cpp index 1a3e3ae1ed72..ada81281dcc6 100644 --- a/SU2_MSH/src/SU2_MSH.cpp +++ b/SU2_MSH/src/SU2_MSH.cpp @@ -285,39 +285,6 @@ int main(int argc, char *argv[]) { grid_adaptation->SetRestart_AdjSolution(config_container[ZONE_0], geo_adapt, config_container[ZONE_0]->GetRestart_AdjFileName()); } - else { - - if (config_container[ZONE_0]->GetKind_Adaptation() == PERIODIC) { - - cout << endl <<"-------------------- Setting the periodic boundaries --------------------" << endl; - - /*--- Set periodic boundary conditions ---*/ - - geometry_container[ZONE_0]->SetPeriodicBoundary(config_container[ZONE_0]); - - /*--- Original grid for debugging purposes ---*/ - - strcpy (file_name, "periodic_original.dat"); geometry_container[ZONE_0]->SetTecPlot(file_name, true); - - /*--- Create a new grid with the right periodic boundary ---*/ - - CGeometry *periodic; periodic = new CPeriodicGeometry(geometry_container[ZONE_0], config_container[ZONE_0]); - periodic->SetPeriodicBoundary(geometry_container[ZONE_0], config_container[ZONE_0]); - periodic->SetMeshFile(geometry_container[ZONE_0], config_container[ZONE_0], config_container[ZONE_0]->GetMesh_Out_FileName()); - - /*--- Output of the grid for debuging purposes ---*/ - - strcpy (file_name, "periodic_halo.dat"); periodic->SetTecPlot(file_name, true); - - } - - if (config_container[ZONE_0]->GetKind_Adaptation() == NONE) { - strcpy (file_name, "original_grid.dat"); - geometry_container[ZONE_0]->SetTecPlot(file_name, true); - geometry_container[ZONE_0]->SetMeshFile(config_container[ZONE_0], config_container[ZONE_0]->GetMesh_Out_FileName()); - } - - } if (rank == MASTER_NODE) cout << endl <<"------------------------- Solver Postprocessing -------------------------" << endl;