Skip to content

Commit

Permalink
Merge pull request #115 from arcaneframework/dev/mab/3D
Browse files Browse the repository at this point in the history 
3D simulations
  • Loading branch information
@mohd-afeef-badri
mohd-afeef-badri authored Feb 3, 2024
2 parents 8458126 + 9953abb commit 4894bc7
Show file tree
Hide file tree
Showing 13 changed files with 1,568 additions and 45 deletions.
4 changes: 4 additions & 0 deletions laplace/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,12 @@ arcane_add_arcane_libraries_to_target(Laplace)
target_include_directories(Laplace PUBLIC . ${CMAKE_CURRENT_BINARY_DIR})
configure_file(Laplace.config ${CMAKE_CURRENT_BINARY_DIR} COPYONLY)
configure_file(Test.laplace.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONLY)
configure_file(Test.laplace.3D.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONLY)
configure_file(Test.laplace.PointDirichlet.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONLY)
configure_file(Test.laplace.PointDirichlet.10K.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONLY)
configure_file(${MSH_DIR}/ring.msh ${CMAKE_CURRENT_BINARY_DIR} COPYONLY)
configure_file(${MSH_DIR}/plancher.msh ${CMAKE_CURRENT_BINARY_DIR} COPYONLY)
configure_file(${MSH_DIR}/L-shape-3D.msh ${CMAKE_CURRENT_BINARY_DIR} COPYONLY)

target_link_libraries(Laplace PUBLIC FemUtils)

Expand Down Expand Up @@ -40,11 +42,13 @@ if(FEMUTILS_HAS_SOLVER_BACKEND_PETSC)
add_test(NAME [laplace]laplace COMMAND Laplace Test.laplace.arc)
endif()

add_test(NAME [laplace]laplace_3D_Dirichlet COMMAND Laplace Test.laplace.3D.arc)
add_test(NAME [laplace]laplace_pointDirichlet COMMAND Laplace Test.laplace.PointDirichlet.arc)

# If parallel part is available, add some tests
if(FEMUTILS_HAS_PARALLEL_SOLVER AND MPIEXEC_EXECUTABLE)
add_test(NAME [laplace]laplace_4pe COMMAND ${MPIEXEC_EXECUTABLE} -n 4 ./Laplace Test.laplace.arc)
add_test(NAME [laplace]laplace_3D_4pe COMMAND ${MPIEXEC_EXECUTABLE} -n 4 ./Laplace Test.laplace.3D.arc)
if(FEMTEST_HAS_GMSH_TEST)
add_test(NAME [laplace]laplace_pointDirichlet_10k_4pe COMMAND ${MPIEXEC_EXECUTABLE} -n 4 ./Laplace Test.laplace.PointDirichlet.10K.arc)
endif()
Expand Down
149 changes: 134 additions & 15 deletions laplace/FemModule.cc
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
//-----------------------------------------------------------------------------
// Copyright 2000-2023 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
// Copyright 2000-2024 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: Apache-2.0
//-----------------------------------------------------------------------------
Expand Down Expand Up @@ -80,14 +80,17 @@ class FemModule
void _updateBoundayConditions();
void _assembleBilinearOperatorTRIA3();
void _assembleBilinearOperatorQUAD4();
void _assembleBilinearOperatorTETRA4();
void _solve();
void _initBoundaryconditions();
void _assembleLinearOperator();
void _applyDirichletBoundaryConditions();
void _checkResultFile();
FixedMatrix<3, 3> _computeElementMatrixTRIA3(Cell cell);
FixedMatrix<4, 4> _computeElementMatrixTETRA4(Cell cell);
FixedMatrix<4, 4> _computeElementMatrixQUAD4(Cell cell);
Real _computeAreaTriangle3(Cell cell);
Real _computeAreaTetra4(Cell cell);
Real _computeAreaQuad4(Cell cell);
Real _computeEdgeLength2(Face face);
Real2 _computeEdgeNormal2(Face face);
Expand Down Expand Up @@ -163,6 +166,8 @@ _doStationarySolve()
// Assemble the FEM bilinear operator (LHS - matrix A)
if (options()->meshType == "QUAD4")
_assembleBilinearOperatorQUAD4();
else if (options()->meshType == "TETRA4")
_assembleBilinearOperatorTETRA4();
else
_assembleBilinearOperatorTRIA3();

Expand All @@ -185,7 +190,7 @@ _getMaterialParameters()
info() << "Get material parameters...";
ElementNodes = 3.;

if (options()->meshType == "QUAD4")
if (options()->meshType == "QUAD4" || options()->meshType == "TETRA4")
ElementNodes = 4.;
}

Expand All @@ -197,7 +202,6 @@ _initBoundaryconditions()
{
info() << "Init boundary conditions...";

info() << "Apply boundary conditions";
_applyDirichletBoundaryConditions();
}

Expand Down Expand Up @@ -442,6 +446,7 @@ _assembleLinearOperator()
}

}

}

/*---------------------------------------------------------------------------*/
Expand All @@ -461,6 +466,26 @@ _computeAreaQuad4(Cell cell)
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

Real FemModule::
_computeAreaTetra4(Cell cell)
{
Real3 m0 = m_node_coord[cell.nodeId(0)];
Real3 m1 = m_node_coord[cell.nodeId(1)];
Real3 m2 = m_node_coord[cell.nodeId(2)];
Real3 m3 = m_node_coord[cell.nodeId(3)];

// Calculate vectors representing edges of the tetrahedron
Real3 v0 = m1 - m0;
Real3 v1 = m2 - m0;
Real3 v2 = m3 - m0;

// Compute volume using scalar triple product
return std::abs(Arcane::math::dot(v0, Arcane::math::cross(v1, v2))) / 6.0;
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

Real FemModule::
_computeAreaTriangle3(Cell cell)
{
Expand Down Expand Up @@ -501,10 +526,76 @@ _computeEdgeNormal2(Face face)
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

FixedMatrix<4, 4> FemModule::
_computeElementMatrixTETRA4(Cell cell)
{
// Get coordinates of the triangle element TETRA4
//------------------------------------------------
// 3 o
// /|\
// / | \
// / | \
// / o 2 \
// / . . \
// o-----------o
// 0 1
//------------------------------------------------
Real3 m0 = m_node_coord[cell.nodeId(0)];
Real3 m1 = m_node_coord[cell.nodeId(1)];
Real3 m2 = m_node_coord[cell.nodeId(2)];
Real3 m3 = m_node_coord[cell.nodeId(3)];

Real volume = _computeAreaTetra4(cell);

// Compute gradients of shape functions
Real3 dPhi0 = Arcane::math::cross(m2 - m1, m1 - m3) ;
Real3 dPhi1 = Arcane::math::cross(m3 - m0, m0 - m2) ;
Real3 dPhi2 = Arcane::math::cross(m1 - m0, m0 - m3) ;
Real3 dPhi3 = Arcane::math::cross(m0 - m1, m1 - m2) ;

// Construct the B-matrix
FixedMatrix<3, 4> b_matrix;
b_matrix(0, 0) = dPhi0.x;
b_matrix(1, 0) = dPhi0.y;
b_matrix(2, 0) = dPhi0.z;

b_matrix(0, 1) = dPhi1.x;
b_matrix(1, 1) = dPhi1.y;
b_matrix(2, 1) = dPhi1.z;

b_matrix(0, 2) = dPhi2.x;
b_matrix(1, 2) = dPhi2.y;
b_matrix(2, 2) = dPhi2.z;

b_matrix(0, 3) = dPhi3.x;
b_matrix(1, 3) = dPhi3.y;
b_matrix(2, 3) = dPhi3.z;

b_matrix.multInPlace(1.0 / (6.0 * volume));

// Compute the element matrix
FixedMatrix<4, 4> int_cdPi_dPj = matrixMultiplication(matrixTranspose(b_matrix), b_matrix);
int_cdPi_dPj.multInPlace(volume);

/*
cout << " Ae \n"
<< "\t" << int_cdPi_dPj(0,0)<<"\t"<< int_cdPi_dPj(0,1)<<"\t"<< int_cdPi_dPj(0,2)<<"\t"<< int_cdPi_dPj(0,3)<<"\n"
<< "\t" << int_cdPi_dPj(1,0)<<"\t"<< int_cdPi_dPj(1,1)<<"\t"<< int_cdPi_dPj(1,2)<<"\t"<< int_cdPi_dPj(1,3)<<"\n"
<< "\t" << int_cdPi_dPj(2,0)<<"\t"<< int_cdPi_dPj(2,1)<<"\t"<< int_cdPi_dPj(2,2)<<"\t"<< int_cdPi_dPj(2,3)<<"\n"
<< "\t" << int_cdPi_dPj(3,0)<<"\t"<< int_cdPi_dPj(3,1)<<"\t"<< int_cdPi_dPj(3,2)<<"\t"<< int_cdPi_dPj(3,3)<<"\n"
<< endl;
*/

return int_cdPi_dPj;
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

FixedMatrix<3, 3> FemModule::
_computeElementMatrixTRIA3(Cell cell)
{
// Get coordiantes of the triangle element TRI3
// Get coordinates of the triangle element TRI3
//------------------------------------------------
// 0 o
// . .
Expand Down Expand Up @@ -536,11 +627,6 @@ _computeElementMatrixTRIA3(Cell cell)
FixedMatrix<3, 3> int_cdPi_dPj = matrixMultiplication(matrixTranspose(b_matrix), b_matrix);
int_cdPi_dPj.multInPlace(area);

//info() << "Cell=" << cell.localId();
//std::cout << " int_cdPi_dPj=";
//int_cdPi_dPj.dump(std::cout);
//std::cout << "\n";

return int_cdPi_dPj;
}

Expand All @@ -550,7 +636,7 @@ _computeElementMatrixTRIA3(Cell cell)
FixedMatrix<4, 4> FemModule::
_computeElementMatrixQUAD4(Cell cell)
{
// Get coordiantes of the quadrangular element QUAD4
// Get coordinates of the quadrangular element QUAD4
//------------------------------------------------
// 1 o . . . . o 0
// . .
Expand Down Expand Up @@ -586,11 +672,6 @@ _computeElementMatrixQUAD4(Cell cell)
FixedMatrix<4, 4> int_cdPi_dPj = matrixMultiplication(matrixTranspose(b_matrix), b_matrix);
int_cdPi_dPj.multInPlace(area);

//info() << "Cell=" << cell.localId();
//std::cout << " int_cdPi_dPj=";
//int_cdPi_dPj.dump(std::cout);
//std::cout << "\n";

return int_cdPi_dPj;
}

Expand Down Expand Up @@ -671,6 +752,44 @@ _assembleBilinearOperatorTRIA3()
}
}
}
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

void FemModule::
_assembleBilinearOperatorTETRA4()
{
auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());

ENUMERATE_ (Cell, icell, allCells()) {
Cell cell = *icell;

auto K_e = _computeElementMatrixTETRA4(cell); // element stiffness matrix

// # assemble elementary matrix into the global one
// # elementary terms are positioned into K according
// # to the rank of associated node in the mesh.nodes list
// for node1 in elem.nodes:
// inode1=elem.nodes.index(node1) # get position of node1 in nodes list
// for node2 in elem.nodes:
// inode2=elem.nodes.index(node2)
// K[node1.rank,node2.rank]=K[node1.rank,node2.rank]+K_e[inode1,inode2]
Int32 n1_index = 0;
for (Node node1 : cell.nodes()) {
Int32 n2_index = 0;
for (Node node2 : cell.nodes()) {
// K[node1.rank,node2.rank]=K[node1.rank,node2.rank]+K_e[inode1,inode2]
Real v = K_e(n1_index, n2_index);
// m_k_matrix(node1.localId(), node2.localId()) += v;
if (node1.isOwn()) {
m_linear_system.matrixAddValue(node_dof.dofId(node1, 0), node_dof.dofId(node2, 0), v);
}
++n2_index;
}
++n1_index;
}
}

}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
Expand Down
33 changes: 33 additions & 0 deletions laplace/Test.laplace.3D.arc
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
<?xml version="1.0"?>
<case codename="Laplace" xml:lang="en" codeversion="1.0">
<arcane>
<title>Sample</title>
<timeloop>LaplaceLoop</timeloop>
</arcane>

<arcane-post-processing>
<output-period>1</output-period>
<output>
<variable>U</variable>
</output>
</arcane-post-processing>

<meshes>
<mesh>
<filename>L-shape-3D.msh</filename>
</mesh>
</meshes>

<fem>
<result-file>test_3D_L-shape.txt</result-file>
<mesh-type>TETRA4</mesh-type>
<dirichlet-boundary-condition>
<surface>bot</surface>
<value>50.0</value>
</dirichlet-boundary-condition>
<dirichlet-boundary-condition>
<surface>bc</surface>
<value>10.0</value>
</dirichlet-boundary-condition>
</fem>
</case>
Loading

0 comments on commit 4894bc7

Please sign in to comment.