Skip to content

Commit

Permalink
test routine for InterfaceDGElasticityIntegrator.
Browse files Browse the repository at this point in the history
  • Loading branch information
dreamer2368 committed Feb 9, 2024
1 parent a032279 commit f8105c4
Show file tree
Hide file tree
Showing 3 changed files with 144 additions and 0 deletions.
1 change: 1 addition & 0 deletions test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ file(COPY inputs/test.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
file(COPY inputs/stokes.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
file(COPY inputs/steady_ns.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
file(COPY meshes/test.2x2.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes)
file(COPY meshes/test.2x1.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes)

file(COPY inputs/test.component.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
file(COPY inputs/stokes.component.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
Expand Down
70 changes: 70 additions & 0 deletions test/dg_integ_mms.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@

#include<gtest/gtest.h>
#include "mms_suite.hpp"
#include "topology_handler.hpp"

using namespace std;
using namespace mfem;
Expand Down Expand Up @@ -34,6 +35,75 @@ TEST(DG_BDR_NORMAL_LF_Test, Test_Tri)
return;
}

// TODO(axel) : test for InterfaceDGElasticityIntegrator.
/* In essence, check the consistency with DGElasticityIntegrator*/
TEST(InterfaceDGElasticityIntegrator, Test_Quad)
{
config = InputParser("inputs/dd_mms.yml");
/* set your own parameters */
const int order = 1;
config.dict_["discretization"]["order"] = order;
config.dict_["mesh"]["filename"] = "meshes/test.2x1.mesh";

SubMeshTopologyHandler *submesh = new SubMeshTopologyHandler;
Mesh *pmesh = submesh->GetGlobalMesh();
const int dim = pmesh->Dimension();
const int udim = dim;
Array<Mesh *> meshes;
TopologyData topol_data;
submesh->ExportInfo(meshes, topol_data);

FiniteElementCollection *fec = new DG_FECollection(order, dim, BasisType::GaussLobatto);
FiniteElementSpace *pfes = new FiniteElementSpace(pmesh, fec, udim);
Array<FiniteElementSpace *> fes(meshes.Size());
for (int k = 0; k < meshes.Size(); k++)
fes[k] = new FiniteElementSpace(meshes[k], fec, udim);

/* will have to initialize this coefficient */
PWConstCoefficient *lambda_cs, *mu_cs;
double alpha, kappa;

/* assemble standard DGElasticityIntegrator Assemble */
BilinearForm pform(pfes);
pform.AddInteriorFaceIntegrator(new DGElasticityIntegrator(*lambda_cs, *mu_cs, alpha, kappa));
pform.Assemble();
pform.Finalize();

/* assemble InterfaceDGElasticityIntegrator Assemble */
InterfaceForm a_itf(meshes, fes, submesh);
a_itf.AddIntefaceIntegrator(new InterfaceDGElasticityIntegrator(alpha, kappa));

Array2D<SparseMatrix *> mats;
mats.SetSize(topol_data.numSub, topol_data.numSub);
for (int i = 0; i < topol_data.numSub; i++)
for (int j = 0; j < topol_data.numSub; j++)
mats(i, j) = new SparseMatrix(fes[i]->GetTrueVSize(), fes[j]->GetTrueVSize());

a_itf.AssembleInterfaceMatrixes(mats);

for (int i = 0; i < topol_data.numSub; i++)
for (int j = 0; j < topol_data.numSub; j++)
mats(i, j)->Finalize();

/*
Assembled matrix from two integrators should have identical non-zero entries.
(but their indices might be different)
Right now their is only one interface with two elements.
Easiest way might be visualizing matrices after changing them into DenseMatrix.
*/
DenseMatrix *pmat = pform.SpMat().ToDenseMatrix();
Array2D<DenseMatrix *> smats(topol_data.numSub, topol_data.numSub);
for (int i = 0; i < topol_data.numSub; i++)
for (int j = 0; j < topol_data.numSub; j++)
smats(i, j) = mats(i, j)->ToDenseMatrix();
// print out or compare the entries...


/* clean up variables */
// delete ..
return;
}

int main(int argc, char* argv[])
{
MPI_Init(&argc, &argv);
Expand Down
73 changes: 73 additions & 0 deletions test/meshes/test.2x1.mesh
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
MFEM mesh v1.0

#
# MFEM Geometry Types (see mesh/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
# TRIANGLE = 2
# SQUARE = 3
# TETRAHEDRON = 4
# CUBE = 5
# PRISM = 6
# PYRAMID = 7
#

dimension
2

elements
2
1 3 0 1 4 3
2 3 1 2 5 4

boundary
6
1 1 0 1
1 1 1 2
3 1 4 3
3 1 5 4
4 1 3 0
2 1 2 5

vertices
6

nodes
FiniteElementSpace
FiniteElementCollection: L2_T1_2D_P3
VDim: 2
Ordering: 1

0 0
0.2763932 0
0.7236068 0
1 0
0 0.2763932
0.2763932 0.2763932
0.7236068 0.2763932
1 0.2763932
0 0.7236068
0.2763932 0.7236068
0.7236068 0.7236068
1 0.7236068
0 1
0.2763932 1
0.7236068 1
1 1
1 0
1.2763932 0
1.7236068 0
2 0
1 0.2763932
1.2763932 0.2763932
1.7236068 0.2763932
2 0.2763932
1 0.7236068
1.2763932 0.7236068
1.7236068 0.7236068
2 0.7236068
1 1
1.2763932 1
1.7236068 1
2 1

0 comments on commit f8105c4

Please sign in to comment.