diff --git a/Common/src/grid_movement_structure.cpp b/Common/src/grid_movement_structure.cpp index 19214653dfec..ffcb71f1d2b4 100644 --- a/Common/src/grid_movement_structure.cpp +++ b/Common/src/grid_movement_structure.cpp @@ -9315,14 +9315,43 @@ void CElasticityMovement::UpdateMultiGrid(CGeometry **geometry, CConfig *config) void CElasticityMovement::SetBoundaryDisplacements(CGeometry *geometry, CConfig *config){ - unsigned short iMarker; - /*--- Get the SU2 module. SU2_CFD will use this routine for dynamically deforming meshes (MARKER_FSI_INTERFACE). ---*/ unsigned short Kind_SU2 = config->GetKind_SU2(); - /*--- First of all, move the FSI interfaces. ---*/ + /*--- Share halo vertex displacements across ranks using the solution vector. + The transfer routines do not do this, i.e. halo vertices have wrong values. ---*/ + + unsigned short iMarker, iDim; + unsigned long iNode, iVertex; + + su2double VarIncrement = 1.0/su2double(config->GetGridDef_Nonlinear_Iter()); + + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + + if ((config->GetMarker_All_ZoneInterface(iMarker) != 0) && (Kind_SU2 == SU2_CFD)) { + + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + + /*--- Get node index ---*/ + iNode = geometry->vertex[iMarker][iVertex]->GetNode(); + + if (geometry->node[iNode]->GetDomain()) { + + /*--- Get the displacement on the vertex ---*/ + for (iDim = 0; iDim < nDim; iDim++) + Solution[iDim] = geometry->vertex[iMarker][iVertex]->GetVarCoord()[iDim] * VarIncrement; + + /*--- Initialize the solution vector ---*/ + LinSysSol.SetBlock(iNode, Solution); + } + } + } + } + StiffMatrix.SendReceive_Solution(LinSysSol, geometry, config); + + /*--- Apply displacement boundary conditions to the FSI interfaces. ---*/ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { if ((config->GetMarker_All_ZoneInterface(iMarker) != 0) && (Kind_SU2 == SU2_CFD)) { @@ -9358,8 +9387,6 @@ void CElasticityMovement::SetClamped_Boundary(CGeometry *geometry, CConfig *conf unsigned long iNode, iVertex; unsigned long iPoint, jPoint; - su2double valJacobian_ij_00 = 0.0; - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Get node index ---*/ @@ -9386,36 +9413,25 @@ void CElasticityMovement::SetClamped_Boundary(CGeometry *geometry, CConfig *conf /*--- Delete the full row for node iNode ---*/ for (jPoint = 0; jPoint < nPoint; jPoint++){ - - /*--- Check whether the block is non-zero ---*/ - valJacobian_ij_00 = StiffMatrix.GetBlock(iNode, jPoint,0,0); - - if (valJacobian_ij_00 != 0.0 ){ - /*--- Set the rest of the row to 0 ---*/ - if (iNode != jPoint) { - StiffMatrix.SetBlock(iNode,jPoint,matrixZeros); - } - /*--- And the diagonal to 1.0 ---*/ - else{ - StiffMatrix.SetBlock(iNode,jPoint,matrixId); - } + if (iNode != jPoint) { + StiffMatrix.SetBlock(iNode,jPoint,matrixZeros); + } + else{ + StiffMatrix.SetBlock(iNode,jPoint,matrixId); } } /*--- Delete the full column for node iNode ---*/ for (iPoint = 0; iPoint < nPoint; iPoint++){ - - /*--- Check whether the block is non-zero ---*/ - valJacobian_ij_00 = StiffMatrix.GetBlock(iPoint, iNode,0,0); - - if (valJacobian_ij_00 != 0.0 ){ - /*--- Set the rest of the row to 0 ---*/ - if (iNode != iPoint) { - StiffMatrix.SetBlock(iPoint,iNode,matrixZeros); - } + if (iNode != iPoint) { + StiffMatrix.SetBlock(iPoint,iNode,matrixZeros); } } } + else { + /*--- Delete the column (iPoint is halo so Send/Recv does the rest) ---*/ + for (iPoint = 0; iPoint < nPoint; iPoint++) StiffMatrix.SetBlock(iPoint,iNode,matrixZeros); + } } } @@ -9424,13 +9440,9 @@ void CElasticityMovement::SetMoving_Boundary(CGeometry *geometry, CConfig *confi unsigned short iDim, jDim; - su2double *VarCoord = NULL; - unsigned long iNode, iVertex; unsigned long iPoint, jPoint; - su2double VarIncrement = 1.0/((su2double)config->GetGridDef_Nonlinear_Iter()); - su2double valJacobian_ij_00 = 0.0; su2double auxJacobian_ij[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; @@ -9440,82 +9452,64 @@ void CElasticityMovement::SetMoving_Boundary(CGeometry *geometry, CConfig *confi iNode = geometry->vertex[val_marker][iVertex]->GetNode(); - /*--- Get the displ;acement on the vertex ---*/ + /*--- Get the displacement of the node ---*/ for (iDim = 0; iDim < nDim; iDim++) - VarCoord = geometry->vertex[val_marker][iVertex]->GetVarCoord(); + Solution[iDim] = LinSysSol.GetBlock(iNode, iDim); if (geometry->node[iNode]->GetDomain()) { - if (nDim == 2) { - Solution[0] = VarCoord[0] * VarIncrement; Solution[1] = VarCoord[1] * VarIncrement; - Residual[0] = VarCoord[0] * VarIncrement; Residual[1] = VarCoord[1] * VarIncrement; - } - else { - Solution[0] = VarCoord[0] * VarIncrement; Solution[1] = VarCoord[1] * VarIncrement; Solution[2] = VarCoord[2] * VarIncrement; - Residual[0] = VarCoord[0] * VarIncrement; Residual[1] = VarCoord[1] * VarIncrement; Residual[2] = VarCoord[2] * VarIncrement; - } - /*--- Initialize the reaction vector ---*/ - LinSysRes.SetBlock(iNode, Residual); - LinSysSol.SetBlock(iNode, Solution); + LinSysRes.SetBlock(iNode, Solution); /*--- STRONG ENFORCEMENT OF THE DISPLACEMENT BOUNDARY CONDITION ---*/ /*--- Delete the full row for node iNode ---*/ for (jPoint = 0; jPoint < nPoint; jPoint++){ - - /*--- Check whether the block is non-zero ---*/ - valJacobian_ij_00 = StiffMatrix.GetBlock(iNode, jPoint,0,0); - - if (valJacobian_ij_00 != 0.0 ){ - /*--- Set the rest of the row to 0 ---*/ - if (iNode != jPoint) { - StiffMatrix.SetBlock(iNode,jPoint,matrixZeros); - } - /*--- And the diagonal to 1.0 ---*/ - else{ - StiffMatrix.SetBlock(iNode,jPoint,matrixId); - } + if (iNode != jPoint) { + StiffMatrix.SetBlock(iNode,jPoint,matrixZeros); + } + else{ + StiffMatrix.SetBlock(iNode,jPoint,matrixId); } } + } - /*--- Delete the columns for a particular node ---*/ - - for (iPoint = 0; iPoint < nPoint; iPoint++){ + /*--- Always delete the iNode column, even for halos ---*/ + for (iPoint = 0; iPoint < nPoint; iPoint++){ - /*--- Check if the term K(iPoint, iNode) is 0 ---*/ - valJacobian_ij_00 = StiffMatrix.GetBlock(iPoint,iNode,0,0); + /*--- Check if the term K(iPoint, iNode) is 0 ---*/ + valJacobian_ij_00 = StiffMatrix.GetBlock(iPoint,iNode,0,0); - /*--- If the node iNode has a crossed dependency with the point iPoint ---*/ - if (valJacobian_ij_00 != 0.0 ){ + /*--- If the node iNode has a crossed dependency with the point iPoint ---*/ + if (valJacobian_ij_00 != 0.0 ){ - /*--- Retrieve the Jacobian term ---*/ - for (iDim = 0; iDim < nDim; iDim++){ - for (jDim = 0; jDim < nDim; jDim++){ - auxJacobian_ij[iDim][jDim] = StiffMatrix.GetBlock(iPoint,iNode,iDim,jDim); - } + /*--- Retrieve the Jacobian term ---*/ + for (iDim = 0; iDim < nDim; iDim++){ + for (jDim = 0; jDim < nDim; jDim++){ + auxJacobian_ij[iDim][jDim] = StiffMatrix.GetBlock(iPoint,iNode,iDim,jDim); } + } - /*--- Multiply by the imposed displacement ---*/ - for (iDim = 0; iDim < nDim; iDim++){ - Residual[iDim] = 0.0; - for (jDim = 0; jDim < nDim; jDim++){ - Residual[iDim] += auxJacobian_ij[iDim][jDim] * VarCoord[jDim]; - } + /*--- Multiply by the imposed displacement ---*/ + for (iDim = 0; iDim < nDim; iDim++){ + Residual[iDim] = 0.0; + for (jDim = 0; jDim < nDim; jDim++){ + Residual[iDim] += auxJacobian_ij[iDim][jDim] * Solution[jDim]; } + } - /*--- For the whole column, except the diagonal term ---*/ - if (iNode != iPoint) { - /*--- The term is substracted from the residual (right hand side) ---*/ - LinSysRes.SubtractBlock(iPoint, Residual); - /*--- The Jacobian term is now set to 0 ---*/ - StiffMatrix.SetBlock(iPoint,iNode,matrixZeros); - } + /*--- For the whole column, except the diagonal term ---*/ + if (iNode != iPoint) { + /*--- The term is substracted from the residual (right hand side) ---*/ + LinSysRes.SubtractBlock(iPoint, Residual); + /*--- The Jacobian term is now set to 0 ---*/ + StiffMatrix.SetBlock(iPoint,iNode,matrixZeros); } } } + } } diff --git a/SU2_CFD/src/solver_direct_elasticity.cpp b/SU2_CFD/src/solver_direct_elasticity.cpp index 9d71a917e692..d68b99a6cbf3 100644 --- a/SU2_CFD/src/solver_direct_elasticity.cpp +++ b/SU2_CFD/src/solver_direct_elasticity.cpp @@ -2513,6 +2513,12 @@ void CFEASolver::BC_Clamped(CGeometry *geometry, CSolver **solver_container, CNu } + } else { + + /*--- Delete the column (iPoint is halo so Send/Recv does the rest) ---*/ + + for (iVar = 0; iVar < nPoint; iVar++) Jacobian.SetBlock(iVar,iPoint,mZeros_Aux); + } } @@ -2647,52 +2653,46 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, CSolver **solver_container, CNu /*--- Delete the full row for node iNode ---*/ for (jPoint = 0; jPoint < nPoint; jPoint++){ - - /*--- Check whether the block is non-zero ---*/ - valJacobian_ij_00 = Jacobian.GetBlock(iNode, jPoint,0,0); - - if (valJacobian_ij_00 != 0.0 ){ - if (iNode != jPoint) { - Jacobian.SetBlock(iNode,jPoint,mZeros_Aux); - } - else{ - Jacobian.SetBlock(iNode,jPoint,mId_Aux); - } + if (iNode != jPoint) { + Jacobian.SetBlock(iNode,jPoint,mZeros_Aux); + } + else{ + Jacobian.SetBlock(iNode,jPoint,mId_Aux); } } - /*--- Delete the columns for a particular node ---*/ + } - for (iPoint = 0; iPoint < nPoint; iPoint++){ + /*--- Always delete the iNode column, even for halos ---*/ - /*--- Check if the term K(iPoint, iNode) is 0 ---*/ - valJacobian_ij_00 = Jacobian.GetBlock(iPoint,iNode,0,0); + for (iPoint = 0; iPoint < nPoint; iPoint++) { - /*--- If the node iNode has a crossed dependency with the point iPoint ---*/ - if (valJacobian_ij_00 != 0.0 ){ + /*--- Check if the term K(iPoint, iNode) is 0 ---*/ + valJacobian_ij_00 = Jacobian.GetBlock(iPoint,iNode,0,0); - /*--- Retrieve the Jacobian term ---*/ - for (iDim = 0; iDim < nDim; iDim++){ - for (jDim = 0; jDim < nDim; jDim++){ - auxJacobian_ij[iDim][jDim] = Jacobian.GetBlock(iPoint,iNode,iDim,jDim); - } - } + /*--- If the node iNode has a crossed dependency with the point iPoint ---*/ + if (valJacobian_ij_00 != 0.0 ){ - /*--- Multiply by the imposed displacement ---*/ - for (iDim = 0; iDim < nDim; iDim++){ - Residual[iDim] = 0.0; - for (jDim = 0; jDim < nDim; jDim++){ - Residual[iDim] += auxJacobian_ij[iDim][jDim] * Disp_Dir[jDim]; - } + /*--- Retrieve the Jacobian term ---*/ + for (iDim = 0; iDim < nDim; iDim++){ + for (jDim = 0; jDim < nDim; jDim++){ + auxJacobian_ij[iDim][jDim] = Jacobian.GetBlock(iPoint,iNode,iDim,jDim); } + } - if (iNode != iPoint) { - /*--- The term is substracted from the residual (right hand side) ---*/ - LinSysRes.SubtractBlock(iPoint, Residual); - /*--- The Jacobian term is now set to 0 ---*/ - Jacobian.SetBlock(iPoint,iNode,mZeros_Aux); + /*--- Multiply by the imposed displacement ---*/ + for (iDim = 0; iDim < nDim; iDim++){ + Residual[iDim] = 0.0; + for (jDim = 0; jDim < nDim; jDim++){ + Residual[iDim] += auxJacobian_ij[iDim][jDim] * Disp_Dir[jDim]; } + } + if (iNode != iPoint) { + /*--- The term is substracted from the residual (right hand side) ---*/ + LinSysRes.SubtractBlock(iPoint, Residual); + /*--- The Jacobian term is now set to 0 ---*/ + Jacobian.SetBlock(iPoint,iNode,mZeros_Aux); } } @@ -2701,7 +2701,6 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, CSolver **solver_container, CNu } - } void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, CNumerics **numerics,