Skip to content

Commit

Permalink
Merge branch 'develop' into tpl/bramwell/update_axom
Browse files Browse the repository at this point in the history
  • Loading branch information
ebchin authored Nov 3, 2023
2 parents aadd698 + 07a71df commit 9363ca9
Show file tree
Hide file tree
Showing 41 changed files with 1,201 additions and 821 deletions.
2 changes: 2 additions & 0 deletions RELEASE-NOTES.md
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@ Changelog](http://keepachangelog.com/en/1.0.0/).
- Renamed getMfemSparseMatrix() to getJacobianSparseMatrix() and getCSRMatrix()
to getJacobianCSRMatrix() to avoid confusion with the separate new MFEM
interface.
- Logging refactor using SLIC macros. Lots of warnings were demoted to DEBUG level.
- Changed various computational geometry routines to return a FaceGeomError enum error handling

### Fixed
- Allow null velocity and response pointers for various use cases
Expand Down
12 changes: 12 additions & 0 deletions src/tests/tribol_common_plane_gap_rate.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -341,6 +341,8 @@ TEST_F( CommonPlaneTest, constant_rate_penetration )
checkPressures( couplingScheme, pressure, 1.E-8, "rate" );
checkForceSense( couplingScheme ); // note: the kinematic and rate contributions are not separated

tribol::finalize();

} // end test 'rate_penetration'

TEST_F( CommonPlaneTest, constant_rate_separation )
Expand Down Expand Up @@ -436,6 +438,8 @@ TEST_F( CommonPlaneTest, constant_rate_separation )
checkPressures( couplingScheme, pressure, 1.E-8, "rate" );
checkForceSense( couplingScheme ); // note: the kinematic and rate contributions aren't separated

tribol::finalize();

} // end test 'rate_separation'

TEST_F( CommonPlaneTest, no_gap_constant_rate_penetration )
Expand Down Expand Up @@ -527,6 +531,8 @@ TEST_F( CommonPlaneTest, no_gap_constant_rate_penetration )
real rate_gap = velZ2 - velZ1;
real pressure = (gap < 0. && rate_gap < 0.) ? penalty * rate_gap : 0.;
checkPressures( couplingScheme, pressure, 1.E-8, "rate" );

tribol::finalize();
}

TEST_F( CommonPlaneTest, percent_rate_penetration )
Expand Down Expand Up @@ -621,6 +627,8 @@ TEST_F( CommonPlaneTest, percent_rate_penetration )
checkPressures( couplingScheme, pressure, 1.E-8, "rate" );
checkForceSense( couplingScheme ); // note: the kinematic and rate contributions are not separated

tribol::finalize();

} // end test 'rate_penetration'

TEST_F( CommonPlaneTest, percent_rate_separation )
Expand Down Expand Up @@ -717,6 +725,8 @@ TEST_F( CommonPlaneTest, percent_rate_separation )
checkPressures( couplingScheme, pressure, 1.E-8, "rate" );
checkForceSense( couplingScheme ); // note: the kinematic and rate contributions aren't separated

tribol::finalize();

} // end test 'rate_separation'

TEST_F( CommonPlaneTest, no_gap_percent_rate_penetration )
Expand Down Expand Up @@ -808,6 +818,8 @@ TEST_F( CommonPlaneTest, no_gap_percent_rate_penetration )
real rate_gap = velZ2 - velZ1;
real pressure = (gap < 0. && rate_gap < 0.) ? penalty * rate_gap : 0.;
checkPressures( couplingScheme, pressure, 1.E-8, "rate" );

tribol::finalize();
}

int main(int argc, char* argv[])
Expand Down
9 changes: 6 additions & 3 deletions src/tests/tribol_common_plane_penalty.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -292,9 +292,6 @@ TEST_F( CommonPlaneTest, penetration_gap_check )
parameters.penalty_ratio = false;
parameters.const_penalty = 1.0;

std::cout << "penetration gap check before slic test before tribolSetupAndUpdate()" << std::endl;
SLIC_INFO("slic test in penetration gap check before tribolSetupAndUpdate()");

int test_mesh_update_err =
this->m_mesh.tribolSetupAndUpdate( tribol::COMMON_PLANE, tribol::PENALTY,
tribol::FRICTIONLESS, false, parameters );
Expand All @@ -311,6 +308,8 @@ TEST_F( CommonPlaneTest, penetration_gap_check )

compareGaps( couplingScheme, gap, 1.E-8, "kinematic_penetration" );

tribol::finalize();

}

TEST_F( CommonPlaneTest, separation_gap_check )
Expand Down Expand Up @@ -372,6 +371,7 @@ TEST_F( CommonPlaneTest, separation_gap_check )

compareGaps( couplingScheme, gap, 1.E-8, "kinematic_separation" );

tribol::finalize();
}

TEST_F( CommonPlaneTest, constant_penalty_check )
Expand Down Expand Up @@ -438,6 +438,7 @@ TEST_F( CommonPlaneTest, constant_penalty_check )
checkPressures( couplingScheme, pressure, 1.E-8 );
checkForceSense( couplingScheme );

tribol::finalize();
}

TEST_F( CommonPlaneTest, element_penalty_check )
Expand Down Expand Up @@ -533,6 +534,7 @@ TEST_F( CommonPlaneTest, element_penalty_check )
checkPressures( couplingScheme, pressure, 1.E-8 );
checkForceSense( couplingScheme );

tribol::finalize();
}

TEST_F( CommonPlaneTest, tied_contact_check )
Expand Down Expand Up @@ -596,6 +598,7 @@ TEST_F( CommonPlaneTest, tied_contact_check )
checkPressures( couplingScheme, pressure, 1.E-8 );
checkForceSense( couplingScheme, true );

tribol::finalize();
}

int main(int argc, char* argv[])
Expand Down
210 changes: 210 additions & 0 deletions src/tests/tribol_comp_geom.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,8 @@ TEST_F( CompGeomTest, common_plane_check )
couplingSchemeManager.getCoupling( 0 );

EXPECT_EQ( userSpecifiedNumOverlaps, couplingScheme->getNumActivePairs() );

tribol::finalize();
}

TEST_F( CompGeomTest, single_mortar_check )
Expand Down Expand Up @@ -178,6 +180,7 @@ TEST_F( CompGeomTest, single_mortar_check )

EXPECT_EQ( userSpecifiedNumOverlaps, couplingScheme->getNumActivePairs() );

tribol::finalize();
}

TEST_F( CompGeomTest, poly_area_centroid_1 )
Expand Down Expand Up @@ -292,6 +295,213 @@ TEST_F( CompGeomTest, poly_area_centroid_2 )
EXPECT_LE( diff_mag, tol );
}

TEST_F( CompGeomTest, 2d_projections_1 )
{
int dim = 2;
int numVerts = 2;
real xy1[dim*numVerts];
real xy2[dim*numVerts];


// Notice how the face vertices are flipped between register mesh and the segment basis eval!

//registerMesh() face 1 x: 1, 0.75, 0.744548
//registerMesh() face 1 y: 1, 0, 0.0902704
//registerMesh() face 1 x: 2, 0.744548, 0.75
//registerMesh() face 1 y: 2, 0.0902704, 0
//(x0,y0) and (x1,y1): (0.744548, 0.0902704), (0.75, 0).
//(px,py): (0.735935, 0.136655)
//SegmentBasis: phi is 1.51907 not between 0. and 1.
//(x0,y0) and (x1,y1): (0.75, 0), (0.744548, 0.0902704).
//(px,py): (0.735935, 0.136655)
//SegmentBasis: phi is 1.51907 not between 0. and 1.

// TODO update to use these face coordinates from testing
// xy1[0] = 0.75;
// xy1[1] = 0. ;
// xy1[2] = 0.744548;
// xy1[3] = 0.0902704;
//
// xy2[0] = 0.744548;
// xy2[1] = 0.0902704;
// xy2[2] = 0.75;
// xy2[3] = 0.;

// this geometry should be in contact
xy1[0] = 0.75;
xy1[1] = 0.;
xy1[2] = 0.727322;
xy1[3] = 0.183039;

xy2[0] = 0.72705;
xy2[1] = 0.182971;
xy2[2] = 0.75;
xy2[3] = 0.;

// compute face normal
real faceNormal1[dim];
real faceNormal2[dim];

real lambdaX1 = xy1[2]-xy1[0];
real lambdaY1 = xy1[3]-xy1[1];

faceNormal1[0] = lambdaY1;
faceNormal1[1] = -lambdaX1;

real lambdaX2 = xy2[2]-xy2[0];
real lambdaY2 = xy2[3]-xy2[1];

faceNormal2[0] = lambdaY2;
faceNormal2[1] = -lambdaX2;

real cxf1[3] = {0., 0., 0.};
real cxf2[3] = {0., 0., 0.};

tribol::VertexAvgCentroid( xy1, dim, numVerts, cxf1[0], cxf1[1], cxf1[2] );
tribol::VertexAvgCentroid( xy2, dim, numVerts, cxf2[0], cxf2[1], cxf2[2] );

// average the vertex averaged centroids of each face to get a pretty good
// estimate of the common plane centroid
real cx[dim];
cx[0] = 0.5*(cxf1[0] + cxf2[0]);
cx[1] = 0.5*(cxf1[1] + cxf2[1]);

real cxProj1[3] = {0., 0., 0.};
real cxProj2[3] = {0., 0., 0.};

tribol::ProjectPointToSegment( cx[0], cx[1], faceNormal1[0], faceNormal1[1],
cxf1[0], cxf1[1], cxProj1[0], cxProj1[1] );
tribol::ProjectPointToSegment( cx[0], cx[1], faceNormal2[0], faceNormal2[1],
cxf2[0], cxf2[1], cxProj2[0], cxProj2[1] );

real diffx1 = std::abs(cxProj1[0] - 0.738595);
real diffy1 = std::abs(cxProj1[1] - 0.0915028);
real diffx2 = std::abs(cxProj2[0] - 0.738591);
real diffy2 = std::abs(cxProj2[1] - 0.0915022);
EXPECT_LE(diffx1, 1.e-6);
EXPECT_LE(diffy1, 1.e-6);
EXPECT_LE(diffx2, 1.e-6);
EXPECT_LE(diffy2, 1.e-6);

// check with call to tribol::update()
tribol::CommType problem_comm = TRIBOL_COMM_WORLD;
tribol::initialize( 2, problem_comm );

real x1[numVerts];
real y1[numVerts];
real x2[numVerts];
real y2[numVerts];

for (int i=0; i<numVerts; ++i)
{
x1[i] = xy1[i*dim];
y1[i] = xy1[i*dim+1];
x2[i] = xy2[i*dim];
y2[i] = xy2[i*dim+1];
}

tribol::IndexType conn1[2] = {0,1};
tribol::IndexType conn2[2] = {0,1};

tribol::registerMesh( 0, 1, 2, &conn1[0], (int)(tribol::LINEAR_EDGE), &x1[0], &y1[0], nullptr );
tribol::registerMesh( 1, 1, 2, &conn2[0], (int)(tribol::LINEAR_EDGE), &x2[0], &y2[0], nullptr );

real fx1[2] = {0., 0.};
real fy1[2] = {0., 0.};
real fx2[2] = {0., 0.};
real fy2[2] = {0., 0.};

tribol::registerNodalResponse( 0, &fx1[0], &fy1[0], nullptr );
tribol::registerNodalResponse( 1, &fx2[0], &fy2[0], nullptr );

tribol::setKinematicConstantPenalty( 0, 1. );
tribol::setKinematicConstantPenalty( 1, 1. );

tribol::registerCouplingScheme( 0, 0, 1,
tribol::SURFACE_TO_SURFACE,
tribol::AUTO,
tribol::COMMON_PLANE,
tribol::FRICTIONLESS,
tribol::PENALTY,
tribol::BINNING_GRID );

tribol::setPenaltyOptions( 0, tribol::KINEMATIC, tribol::KINEMATIC_CONSTANT );
tribol::setContactPenFrac(0.5);
tribol::setContactAreaFrac(1.e-4);

// TODO check penetration and overlap tolerance with what is being used in host-code

double dt = 1.;
int update_err = tribol::update( 1, 1., dt );

EXPECT_EQ( update_err, 0 );

}

TEST_F( CompGeomTest, 2d_projections_2 )
{
int dim = 2;
int numVerts = 2;
real xy1[dim*numVerts];
real xy2[dim*numVerts];

// face coordinates from testing
xy1[0] = 0.75;
xy1[1] = 0.;
xy1[2] = 0.727322;
xy1[3] = 0.183039;

xy2[0] = 0.727322;
xy2[1] = 0.183039;
xy2[2] = 0.75;
xy2[3] = 0.;

// compute face normal
real faceNormal1[dim];
real faceNormal2[dim];

real lambdaX1 = xy1[2]-xy1[0];
real lambdaY1 = xy1[3]-xy1[1];

faceNormal1[0] = lambdaY1;
faceNormal1[1] = -lambdaX1;

real lambdaX2 = xy2[2]-xy2[0];
real lambdaY2 = xy2[3]-xy2[1];

faceNormal2[0] = lambdaY2;
faceNormal2[1] = -lambdaX2;

real cxf1[3] = {0., 0., 0.};
real cxf2[3] = {0., 0., 0.};

tribol::VertexAvgCentroid( xy1, dim, numVerts, cxf1[0], cxf1[1], cxf1[2] );
tribol::VertexAvgCentroid( xy2, dim, numVerts, cxf2[0], cxf2[1], cxf2[2] );

// average the vertex averaged centroids of each face to get a pretty good
// estimate of the common plane centroid
real cx[dim];
cx[0] = 0.5*(cxf1[0] + cxf2[0]);
cx[1] = 0.5*(cxf1[1] + cxf2[1]);

real cxProj1[3] = {0., 0., 0.};
real cxProj2[3] = {0., 0., 0.};

tribol::ProjectPointToSegment( cx[0], cx[1], faceNormal1[0], faceNormal1[1],
cxf1[0], cxf1[1], cxProj1[0], cxProj1[1] );
tribol::ProjectPointToSegment( cx[0], cx[1], faceNormal2[0], faceNormal2[1],
cxf2[0], cxf2[1], cxProj2[0], cxProj2[1] );

real diffx1 = std::abs(cxProj1[0] - cx[0]);
real diffy1 = std::abs(cxProj1[1] - cx[1]);
real diffx2 = std::abs(cxProj2[0] - cx[0]);
real diffy2 = std::abs(cxProj2[1] - cx[1]);
EXPECT_LE(diffx1, 1.e-6);
EXPECT_LE(diffy1, 1.e-6);
EXPECT_LE(diffx2, 1.e-6);
EXPECT_LE(diffy2, 1.e-6);
}

int main(int argc, char* argv[])
{
int result = 0;
Expand Down
Loading

0 comments on commit 9363ca9

Please sign in to comment.