Skip to content

Fix given displacement BC's of FEA solver and CElasticityMovement #658

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pcarruscag merged 5 commits into from
Mar 1, 2019
Merged

Fix given displacement BC's of FEA solver and CElasticityMovement #658

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0a7a6be
delete columns of halo points on given displacement boundaries
pcarruscag Feb 22, 2019
6eb57be
apply same fix to CElasticityMovement
pcarruscag Feb 22, 2019
2f496ac
fix issue in CElasticityMovement, due to halo vertices not having cor…
pcarruscag Feb 25, 2019
02d9354
revert travis
pcarruscag Feb 26, 2019
3b7cab8
Merge branch 'develop' into fix_fea_disp_bc
rsanfer Mar 1, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
158 changes: 76 additions & 82 deletions Common/src/grid_movement_structure.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)) {
Expand Down Expand Up @@ -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 ---*/
Expand All @@ -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) {
Copy link
Member

@economon economon Feb 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are several places where you do a similar operation as to what is done here.. have you considered reusing something like StiffMatrix.DeleteValsRowi() or creating a new version that suits what you are doing here within the matrix class to save some code by calling it instead?

Copy link
Member Author

@pcarruscag pcarruscag Feb 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes definitely, especially since these operations can be performed more efficiently if the sparsity pattern is know. And there are other areas too, for example we currently add the contribution of the mass matrix to the stiffness matrix in a similar block by block way. But I would prefer to address all of those in one PR and avoid creating conflicts between ongoing PR's.

Copy link
Contributor

@rsanfer rsanfer Mar 1, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think creating a StiffMatrix.DeleteValsColi() would do the job. There is some room for code saving in all of this implementation, but I agree with @pcarruscag, that could be done in an additional PR.

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);
}
}

}
Expand All @@ -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}};

Expand All @@ -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);
}
}
}

}

}
Expand Down
69 changes: 34 additions & 35 deletions SU2_CFD/src/solver_direct_elasticity.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);

}

}
Expand Down Expand Up @@ -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);
}

}
Expand All @@ -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,
Expand Down