diff --git a/laplace/CMakeLists.txt b/laplace/CMakeLists.txt index 742e21e..d746d0e 100644 --- a/laplace/CMakeLists.txt +++ b/laplace/CMakeLists.txt @@ -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) @@ -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() diff --git a/laplace/FemModule.cc b/laplace/FemModule.cc index 1554406..a1af7bf 100644 --- a/laplace/FemModule.cc +++ b/laplace/FemModule.cc @@ -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 //----------------------------------------------------------------------------- @@ -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); @@ -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(); @@ -185,7 +190,7 @@ _getMaterialParameters() info() << "Get material parameters..."; ElementNodes = 3.; - if (options()->meshType == "QUAD4") + if (options()->meshType == "QUAD4" || options()->meshType == "TETRA4") ElementNodes = 4.; } @@ -197,7 +202,6 @@ _initBoundaryconditions() { info() << "Init boundary conditions..."; - info() << "Apply boundary conditions"; _applyDirichletBoundaryConditions(); } @@ -442,6 +446,7 @@ _assembleLinearOperator() } } + } /*---------------------------------------------------------------------------*/ @@ -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) { @@ -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 // . . @@ -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; } @@ -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 // . . @@ -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; } @@ -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; + } + } + +} /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ diff --git a/laplace/Test.laplace.3D.arc b/laplace/Test.laplace.3D.arc new file mode 100644 index 0000000..4e4cf29 --- /dev/null +++ b/laplace/Test.laplace.3D.arc @@ -0,0 +1,33 @@ + + + + Sample + LaplaceLoop + + + + 1 + + U + + + + + + L-shape-3D.msh + + + + + test_3D_L-shape.txt + TETRA4 + + bot + 50.0 + + + bc + 10.0 + + + diff --git a/laplace/tests/test_3D_L-shape.txt b/laplace/tests/test_3D_L-shape.txt new file mode 100644 index 0000000..478234e --- /dev/null +++ b/laplace/tests/test_3D_L-shape.txt @@ -0,0 +1,108 @@ +1 50 +2 50 +3 50 +4 50 +5 34.2593872958221 +6 34.2608949043884 +7 23.8720119291307 +8 23.8873187523708 +9 10 +10 10 +11 34.1601810418892 +12 34.1827229081516 +13 10 +14 10 +15 50 +16 50 +17 50 +18 50 +19 45.2028558172978 +20 40.6333703464156 +21 37.2395745599342 +22 34.9737846264254 +23 35.0072173875521 +24 37.134257114298 +25 40.7820018916361 +26 45.2007928303114 +27 33.5786932590157 +28 31.4756974308377 +29 28.2841911519725 +30 25.4464084903849 +31 25.2264650682491 +32 28.3381535655469 +33 31.4447376517162 +34 33.4498143562111 +35 22.9218550251377 +36 19.0999098513826 +37 19.0872134345656 +38 22.8542563642675 +39 20.9996810836613 +40 27.2722503632292 +41 27.2169472666949 +42 20.964464598427 +43 45.3979446750982 +44 40.2491887324268 +45 40.2795431030281 +46 45.3965792415293 +47 33.0433314942744 +48 24.4872938447329 +49 39.1042042712379 +50 46.3067426294272 +51 35.1792809363659 +52 31.3897328935288 +53 38.2682174266064 +54 19.5326855591507 +55 22.6610017833454 +56 27.0023572384265 +57 42.8556701985193 +58 35.2398936212365 +59 30.697987423999 +60 43.4472961203008 +61 34.1795221790386 +62 24.071786753317 +63 33.105115931429 +64 46.7759289525899 +65 46.8482992161135 +66 28.4511490176226 +67 32.9949839409043 +68 24.3370306692837 +69 39.0599737831777 +70 46.3278937029904 +71 35.2027525699533 +72 31.0599196059805 +73 38.4337394406492 +74 19.636107971884 +75 22.6892618937892 +76 26.7922809991088 +77 42.8609780060907 +78 35.2189019206004 +79 30.675528143092 +80 43.5497442510856 +81 34.2491167224947 +82 24.2415561954037 +83 33.12246040574 +84 46.7796687592443 +85 46.8711232816625 +86 28.3143286597566 +87 34.0377979268469 +88 24.344078119591 +89 30.0497730651866 +90 32.6801964441339 +91 26.7642486298875 +92 47.5708470736884 +93 34.4913830247977 +94 38.8014446550279 +95 42.8821368633878 +96 35.9230090555142 +97 50 +98 50 +99 50 +100 47.2521483735211 +101 42.9432136851237 +102 37.8712709880355 +103 29.5708315612763 +104 24.3692075012447 +105 18.235797144708 +106 15.8508073684188 +107 21.3575688480501 +108 23.4549617696854 diff --git a/meshes/msh/L-shape-3D.geo b/meshes/msh/L-shape-3D.geo new file mode 100644 index 0000000..fae841f --- /dev/null +++ b/meshes/msh/L-shape-3D.geo @@ -0,0 +1,97 @@ +//----------------------------------------------------------------------------- +// +// Name : L-shape-3D.geo +// Author : Mohd Afeef BADRI +// Date : 03 / Feb / 2024 +// +// ---------------------------------------------------------------------------- +// Comment : Simple problem L-shape mesh in 3D. The boundaries of the mesh +// are named "bot", "bc", "other", besides volume is named "vol". +// +// Usage : gmsh L-shape-3D.geo -3 +// +// ---------------------------------------------------------------------------- + +ls=100.; +baselength=250.0; +basewidth=100.0; +baseheight=500.0; +toplength=500.0; +Lheight=250.0; +Lwidth=250.0; + +Point(newp) = {0, 0, 0, ls}; +Point(newp) = {baselength, 0, 0, ls}; +Point(newp) = {0, basewidth, 0, ls}; +Point(newp) = {baselength, basewidth, 0, ls}; +Point(newp) = {0, 0, baseheight, ls}; +Point(newp) = {0, basewidth, baseheight, ls}; +Point(newp) = {toplength, 0, baseheight, ls}; +Point(newp) = {toplength, basewidth, baseheight, ls}; +Point(newp) = {toplength, 0, Lheight, ls}; +Point(newp) = {toplength, basewidth, Lheight, ls}; +Point(newp) = {Lwidth, 0, Lheight, ls}; +Point(newp) = {Lwidth, basewidth, Lheight, ls}; +Point(newp) = {Lwidth+(Lwidth-30.), 0, Lheight, ls}; +Point(newp) = {Lwidth+(Lwidth-30.), basewidth, Lheight, ls}; + + +Line(1) = {1, 3}; +Line(2) = {3, 4}; +Line(3) = {4, 2}; +Line(4) = {2, 1}; +Line(5) = {1, 5}; +Line(6) = {5, 6}; +Line(7) = {6, 3}; +Line(8) = {6, 8}; +Line(9) = {8, 7}; +Line(10) = {7, 5}; +Line(11) = {8, 10}; +Line(12) = {10, 9}; +Line(13) = {9, 7}; +Line(14) = {9, 13}; +Line(15) = {13, 14}; +Line(16) = {14, 10}; +Line(17) = {13, 11}; +Line(18) = {11, 12}; +Line(19) = {12, 14}; +Line(20) = {2, 11}; +Line(21) = {12, 4}; + +Curve Loop(1) = {5, -10, -13, 14, 17, -20, 4}; +Plane Surface(1) = {1}; + +Curve Loop(2) = {8, 11, -16, -19, 21, -2, -7}; +Plane Surface(2) = {2}; + +Curve Loop(3) = {8, 9, 10, 6}; +Plane Surface(3) = {3}; + +Curve Loop(4) = {6, 7, -1, 5}; +Plane Surface(4) = {4}; + +Curve Loop(5) = {2, 3, 4, 1}; +Plane Surface(5) = {5}; + +Curve Loop(6) = {18, 21, 3, 20}; +Plane Surface(6) = {6}; + +Curve Loop(7) = {19, -15, 17, 18}; +Plane Surface(7) = {7}; + +Curve Loop(8) = {12, 14, 15, 16}; +Plane Surface(8) = {8}; + +Curve Loop(9) = {9, -13, -12, -11}; +Plane Surface(9) = {9}; + +Surface Loop(1) = {1, 4, 3, 2, 9, 8, 7, 6, 5}; +Volume(1) = {1}; + +Physical Volume("vol") = {1}; +Physical Surface("bot") = {5}; +Physical Surface("bc") = {8}; +Physical Surface("other") = {1, 7, 9, 3, 2, 6, 4}; + +Mesh.MshFileVersion = 4.1; +Mesh 2; diff --git a/meshes/msh/L-shape-3D.msh b/meshes/msh/L-shape-3D.msh new file mode 100644 index 0000000..9f64db8 --- /dev/null +++ b/meshes/msh/L-shape-3D.msh @@ -0,0 +1,797 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +4 +2 2 "bot" +2 3 "bc" +2 4 "other" +3 1 "vol" +$EndPhysicalNames +$Entities +14 21 9 1 +1 0 0 0 0 +2 250 0 0 0 +3 0 100 0 0 +4 250 100 0 0 +5 0 0 500 0 +6 0 100 500 0 +7 500 0 500 0 +8 500 100 500 0 +9 500 0 250 0 +10 500 100 250 0 +11 250 0 250 0 +12 250 100 250 0 +13 470 0 250 0 +14 470 100 250 0 +1 0 0 0 0 100 0 0 2 1 -3 +2 0 100 0 250 100 0 0 2 3 -4 +3 250 0 0 250 100 0 0 2 4 -2 +4 0 0 0 250 0 0 0 2 2 -1 +5 0 0 0 0 0 500 0 2 1 -5 +6 0 0 500 0 100 500 0 2 5 -6 +7 0 100 0 0 100 500 0 2 6 -3 +8 0 100 500 500 100 500 0 2 6 -8 +9 500 0 500 500 100 500 0 2 8 -7 +10 0 0 500 500 0 500 0 2 7 -5 +11 500 100 250 500 100 500 0 2 8 -10 +12 500 0 250 500 100 250 0 2 10 -9 +13 500 0 250 500 0 500 0 2 9 -7 +14 470 0 250 500 0 250 0 2 9 -13 +15 470 0 250 470 100 250 0 2 13 -14 +16 470 100 250 500 100 250 0 2 14 -10 +17 250 0 250 470 0 250 0 2 13 -11 +18 250 0 250 250 100 250 0 2 11 -12 +19 250 100 250 470 100 250 0 2 12 -14 +20 250 0 0 250 0 250 0 2 2 -11 +21 250 100 0 250 100 250 0 2 12 -4 +1 0 0 0 500 0 500 1 4 7 5 -10 -13 14 17 -20 4 +2 0 100 0 500 100 500 1 4 7 8 11 -16 -19 21 -2 -7 +3 0 0 500 500 100 500 1 4 4 8 9 10 6 +4 0 0 0 0 100 500 1 4 4 6 7 -1 5 +5 0 0 0 250 100 0 1 2 4 2 3 4 1 +6 250 0 0 250 100 250 1 4 4 18 21 3 20 +7 250 0 250 470 100 250 1 4 4 19 -15 17 18 +8 470 0 250 500 100 250 1 3 4 12 14 15 16 +9 500 0 250 500 100 500 1 4 4 9 -13 -12 -11 +1 0 0 0 500 100 500 1 1 9 1 4 3 2 9 8 7 6 5 +$EndEntities +$Nodes +36 108 1 108 +0 1 0 1 +1 +0 0 0 +0 2 0 1 +2 +250 0 0 +0 3 0 1 +3 +0 100 0 +0 4 0 1 +4 +250 100 0 +0 5 0 1 +5 +0 0 500 +0 6 0 1 +6 +0 100 500 +0 7 0 1 +7 +500 0 500 +0 8 0 1 +8 +500 100 500 +0 9 0 1 +9 +500 0 250 +0 10 0 1 +10 +500 100 250 +0 11 0 1 +11 +250 0 250 +0 12 0 1 +12 +250 100 250 +0 13 0 1 +13 +470 0 250 +0 14 0 1 +14 +470 100 250 +1 2 0 2 +15 +16 +83.33333333325572 100 0 +166.6666666666362 100 0 +1 4 0 2 +17 +18 +166.6666666666722 0 0 +83.3333333333278 0 0 +1 5 0 4 +19 +20 +21 +22 +0 0 99.99999999989699 +0 0 199.999999999857 +0 0 299.9999999999268 +0 0 399.9999999999634 +1 7 0 4 +23 +24 +25 +26 +0 100 400.0000000000067 +0 100 300.0000000000133 +0 100 199.9999999999868 +0 100 99.99999999999335 +1 8 0 4 +27 +28 +29 +30 +99.99999999989699 100 500 +199.999999999857 100 500 +299.9999999999268 100 500 +399.9999999999634 100 500 +1 10 0 4 +31 +32 +33 +34 +400.0000000000067 0 500 +300.0000000000133 0 500 +199.9999999999868 0 500 +99.99999999999335 0 500 +1 11 0 2 +35 +36 +500 100 416.6666666666667 +500 100 333.3333333333334 +1 13 0 2 +37 +38 +500 0 333.3333333333333 +500 0 416.6666666666666 +1 17 0 2 +39 +40 +396.6666666666925 0 250 +323.3333333333852 0 250 +1 19 0 2 +41 +42 +323.3333333333852 100 250 +396.6666666666925 100 250 +1 20 0 2 +43 +44 +250 0 83.33333333325572 +250 0 166.6666666666362 +1 21 0 2 +45 +46 +250 100 166.6666666666722 +250 100 83.3333333333278 +2 1 0 20 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +58 +59 +60 +61 +62 +63 +64 +65 +66 +148.5850198415951 0 423.6705016648935 +366.3236461064289 0 310.8358271063244 +175.2688313872357 0 205.8354262330073 +125 0 72.16878364869841 +175.539603548702 0 287.7019508826543 +248.2888919674633 0 334.5455510363292 +92.94808004135396 0 247.6700167963301 +428.8449971857077 0 303.4051727183418 +426.8088049893536 0 365.4585324089549 +336.1941983855479 0 393.2897310891918 +90.11370088011795 0 145.4067790969926 +96.43213200289843 0 348.6920703833276 +236.3025416710615 0 417.7521479811378 +179.8963830505446 0 122.4025371571095 +73.20508075687556 0 426.7949192431129 +432.6006006749816 0 435.0829860329627 +181.0296378063441 0 362.4724443896685 +188.9957660359141 0 61.00423396406664 +59.68940684268915 0 63.51511254911762 +304.8280139585651 0 307.7342218463691 +2 2 0 20 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +81 +82 +83 +84 +85 +86 +148.5850198415466 100 423.6705016649061 +366.3236461064291 100 310.8358271063258 +175.2688313872339 100 205.8354262330171 +124.999999999946 100 72.1687836486828 +175.53960354871 100 287.7019508826634 +248.2888919674608 100 334.5455510363393 +92.94808004135953 100 247.6700167963611 +428.844997185708 100 303.4051727183432 +426.8088049893518 100 365.4585324089589 +336.1941983855348 100 393.2897310891962 +90.11370088010113 100 145.4067790970275 +96.43213200288376 100 348.6920703833689 +236.3025416710133 100 417.7521479811526 +179.8963830505299 100 122.4025371571366 +73.20508075681632 100 426.7949192431578 +432.60060067497 100 435.0829860329644 +181.0296378063229 100 362.4724443896861 +188.9957660358984 100 61.00423396409199 +59.68940684266056 100 63.51511254914074 +304.8280139585619 100 307.7342218463722 +2 3 0 5 +87 +88 +89 +90 +91 +49.99999999994849 49.99999999999999 500 +450.0000000000034 50.00000000000001 500 +249.9999999999459 50 500 +149.999999999877 50.00000000004818 500 +350.0000000000101 49.99999999997839 500 +2 4 0 5 +92 +93 +94 +95 +96 +0 49.99999999999999 49.99999999994849 +0 50.00000000000001 450.0000000000034 +0 50 249.999999999946 +0 49.99999999995183 149.999999999877 +0 50.00000000002159 350.0000000000101 +2 5 0 3 +97 +98 +99 +50.17006802717789 50.00000000000306 0 +124.999999999946 50.00000000003003 0 +199.8299319727824 50.00000000000306 0 +2 6 0 3 +100 +101 +102 +250 49.99999999999694 50.17006802717789 +250 49.99999999996997 124.999999999946 +250 49.99999999999694 199.8299319727824 +2 7 0 3 +103 +104 +105 +301.3333333333618 50 250 +360.000000000039 50 250 +418.6666666666848 50 250 +2 8 0 0 +2 9 0 3 +106 +107 +108 +500.0000000000001 50.00000000000001 300.1700680272108 +500.0000000000001 50.00000000000009 375 +500.0000000000001 50 449.8299319727892 +3 1 0 0 +$EndNodes +$Elements +10 471 1 471 +2 1 2 61 +1 18 1 65 +2 1 19 65 +3 2 17 64 +4 43 2 64 +5 22 5 61 +6 5 34 61 +7 31 7 62 +8 7 38 62 +9 9 13 37 +10 40 11 66 +11 11 44 49 +12 11 49 51 +13 11 51 52 +14 11 52 66 +15 37 13 54 +16 13 39 54 +17 17 18 50 +18 17 50 64 +19 50 18 65 +20 19 20 57 +21 19 57 65 +22 20 21 53 +23 20 53 57 +24 21 22 58 +25 53 21 58 +26 58 22 61 +27 32 31 56 +28 56 31 62 +29 33 32 59 +30 32 56 59 +31 34 33 47 +32 47 33 59 +33 34 47 61 +34 38 37 55 +35 37 54 55 +36 38 55 62 +37 39 40 48 +38 39 48 54 +39 48 40 66 +40 44 43 60 +41 60 43 64 +42 49 44 60 +43 47 58 61 +44 58 47 63 +45 47 59 63 +46 54 48 55 +47 55 48 56 +48 56 48 66 +49 51 49 53 +50 53 49 57 +51 57 49 60 +52 50 57 60 +53 57 50 65 +54 50 60 64 +55 52 51 63 +56 51 53 58 +57 51 58 63 +58 56 52 59 +59 52 56 66 +60 59 52 63 +61 55 56 62 +2 2 2 61 +62 15 3 85 +63 3 26 85 +64 4 16 84 +65 46 4 84 +66 23 6 81 +67 6 27 81 +68 30 8 82 +69 8 35 82 +70 10 14 36 +71 41 12 86 +72 12 45 69 +73 12 69 71 +74 12 71 72 +75 12 72 86 +76 36 14 74 +77 14 42 74 +78 16 15 70 +79 70 15 85 +80 16 70 84 +81 24 23 78 +82 78 23 81 +83 25 24 73 +84 73 24 78 +85 26 25 77 +86 25 73 77 +87 26 77 85 +88 27 28 67 +89 27 67 81 +90 28 29 79 +91 67 28 79 +92 29 30 76 +93 29 76 79 +94 76 30 82 +95 35 36 75 +96 35 75 82 +97 36 74 75 +98 42 41 68 +99 68 41 86 +100 42 68 74 +101 45 46 80 +102 69 45 80 +103 80 46 84 +104 67 78 81 +105 78 67 83 +106 67 79 83 +107 74 68 75 +108 75 68 76 +109 76 68 86 +110 71 69 73 +111 73 69 77 +112 77 69 80 +113 70 77 80 +114 77 70 85 +115 70 80 84 +116 72 71 83 +117 71 73 78 +118 71 78 83 +119 76 72 79 +120 72 76 86 +121 79 72 83 +122 75 76 82 +2 3 2 20 +123 5 6 87 +124 34 5 87 +125 6 27 87 +126 8 7 88 +127 7 31 88 +128 30 8 88 +129 27 28 90 +130 87 27 90 +131 28 29 89 +132 33 28 89 +133 28 33 90 +134 29 30 91 +135 29 32 89 +136 32 29 91 +137 30 88 91 +138 31 32 91 +139 88 31 91 +140 32 33 89 +141 33 34 90 +142 34 87 90 +2 4 2 20 +143 3 1 92 +144 1 19 92 +145 26 3 92 +146 5 6 93 +147 22 5 93 +148 6 23 93 +149 19 20 95 +150 92 19 95 +151 20 21 94 +152 25 20 94 +153 20 25 95 +154 21 22 96 +155 21 24 94 +156 24 21 96 +157 22 93 96 +158 23 24 96 +159 93 23 96 +160 24 25 94 +161 25 26 95 +162 26 92 95 +2 5 2 12 +163 1 3 97 +164 18 1 97 +165 4 2 99 +166 2 17 99 +167 3 15 97 +168 16 4 99 +169 15 16 98 +170 97 15 98 +171 98 16 99 +172 17 18 98 +173 17 98 99 +174 18 97 98 +2 6 2 12 +175 4 2 100 +176 2 43 100 +177 46 4 100 +178 11 12 102 +179 44 11 102 +180 12 45 102 +181 43 44 101 +182 100 43 101 +183 101 44 102 +184 45 46 101 +185 45 101 102 +186 46 100 101 +2 7 2 12 +187 11 12 103 +188 40 11 103 +189 12 41 103 +190 14 13 105 +191 13 39 105 +192 42 14 105 +193 39 40 104 +194 39 104 105 +195 40 103 104 +196 41 42 104 +197 103 41 104 +198 104 42 105 +2 8 2 2 +199 13 10 9 +200 13 14 10 +2 9 2 12 +201 8 7 108 +202 7 38 108 +203 35 8 108 +204 9 10 106 +205 37 9 106 +206 10 36 106 +207 36 35 107 +208 107 35 108 +209 106 36 107 +210 38 37 107 +211 37 106 107 +212 38 107 108 +3 1 4 259 +213 94 53 73 78 +214 25 20 53 94 +215 94 53 78 58 +216 78 73 94 24 +217 53 94 21 58 +218 24 94 78 21 +219 25 53 73 94 +220 21 94 78 58 +221 80 60 101 102 +222 89 59 79 56 +223 55 74 106 54 +224 55 54 48 74 +225 80 70 50 77 +226 50 70 85 77 +227 74 54 48 105 +228 78 71 83 63 +229 60 80 77 69 +230 85 50 98 70 +231 28 33 79 89 +232 65 57 50 85 +233 74 48 55 68 +234 84 98 50 70 +235 63 83 78 47 +236 50 60 80 77 +237 53 57 95 77 +238 101 60 44 102 +239 45 80 101 102 +240 77 95 53 73 +241 50 85 98 65 +242 33 59 79 89 +243 32 29 56 89 +244 63 79 83 47 +245 50 64 98 84 +246 20 25 53 95 +247 24 78 96 21 +248 89 59 56 32 +249 29 89 79 76 +250 56 79 89 76 +251 25 53 95 73 +252 97 65 98 85 +253 78 96 21 58 +254 70 50 84 80 +255 60 80 69 102 +256 95 77 25 73 +257 20 53 57 95 +258 67 78 83 47 +259 11 71 102 51 +260 71 11 102 12 +261 55 75 74 68 +262 75 107 55 106 +263 48 104 68 105 +264 74 55 106 75 +265 99 98 64 84 +266 64 100 84 101 +267 67 79 90 47 +268 64 50 60 84 +269 77 49 60 69 +270 68 74 48 105 +271 71 49 102 51 +272 49 77 60 57 +273 79 67 83 47 +274 57 65 95 85 +275 71 78 51 63 +276 73 49 77 69 +277 73 71 78 51 +278 55 48 56 68 +279 53 77 49 57 +280 77 53 49 73 +281 104 103 66 86 +282 60 50 57 77 +283 29 56 89 76 +284 49 71 102 69 +285 52 83 71 63 +286 66 48 104 68 +287 50 60 84 80 +288 78 53 73 51 +289 106 54 74 105 +290 55 76 75 68 +291 51 11 49 102 +292 71 69 12 102 +293 85 57 50 77 +294 49 60 69 102 +295 60 84 101 64 +296 86 66 104 68 +297 53 78 58 51 +298 78 58 51 63 +299 52 79 83 63 +300 76 66 72 86 +301 47 79 90 59 +302 48 66 56 68 +303 47 63 79 59 +304 63 78 58 47 +305 55 76 82 75 +306 86 76 66 68 +307 67 81 78 47 +308 84 60 101 80 +309 71 73 49 51 +310 72 11 12 103 +311 103 72 66 86 +312 33 28 79 90 +313 85 95 57 77 +314 59 52 79 56 +315 71 51 52 63 +316 29 32 56 91 +317 11 71 51 52 +318 69 73 49 71 +319 72 11 71 12 +320 59 90 47 33 +321 67 79 28 90 +322 72 79 52 56 +323 92 19 65 85 +324 95 65 19 85 +325 62 76 91 82 +326 59 79 90 33 +327 55 107 37 106 +328 36 107 75 106 +329 53 49 73 51 +330 76 62 55 82 +331 11 72 52 103 +332 56 52 72 66 +333 72 83 71 52 +334 81 67 90 47 +335 60 49 44 102 +336 69 80 45 102 +337 61 78 81 47 +338 11 72 71 52 +339 18 65 98 97 +340 15 97 98 85 +341 74 36 75 106 +342 55 37 54 106 +343 63 52 79 59 +344 55 82 107 75 +345 79 83 72 52 +346 91 56 29 76 +347 58 78 61 47 +348 66 52 72 103 +349 64 17 98 99 +350 84 99 98 16 +351 100 84 101 46 +352 101 64 100 43 +353 95 92 85 19 +354 92 95 85 26 +355 76 55 56 68 +356 96 78 81 61 +357 61 81 90 47 +358 72 79 56 76 +359 76 62 91 56 +360 82 55 107 62 +361 107 62 108 82 +362 62 55 56 76 +363 72 56 66 76 +364 104 48 39 105 +365 104 42 68 105 +366 96 78 61 58 +367 18 65 50 98 +368 98 70 15 85 +369 76 56 66 68 +370 90 87 81 34 +371 87 90 81 27 +372 86 104 103 41 +373 103 66 40 104 +374 100 99 64 84 +375 34 87 81 61 +376 93 96 81 61 +377 34 81 90 61 +378 95 85 26 77 +379 57 65 19 95 +380 64 17 50 98 +381 70 84 98 16 +382 12 72 103 86 +383 66 11 52 103 +384 96 81 23 93 +385 96 93 22 61 +386 35 107 108 82 +387 107 62 38 108 +388 39 48 54 105 +389 68 42 74 105 +390 106 74 14 105 +391 84 80 101 46 +392 101 60 64 43 +393 67 27 81 90 +394 90 34 61 47 +395 20 53 94 21 +396 94 73 25 24 +397 13 54 106 105 +398 68 104 86 41 +399 48 40 66 104 +400 88 91 30 62 +401 88 91 62 31 +402 50 17 18 98 +403 98 70 16 15 +404 97 92 65 85 +405 30 88 62 82 +406 30 62 91 82 +407 89 33 59 32 +408 79 28 89 29 +409 14 13 106 105 +410 104 48 40 39 +411 42 104 68 41 +412 78 96 81 23 +413 61 22 96 58 +414 91 76 30 82 +415 62 91 56 31 +416 82 35 107 75 +417 62 55 107 38 +418 21 96 22 58 +419 96 78 24 23 +420 45 101 80 46 +421 101 44 60 43 +422 77 26 95 25 +423 19 20 57 95 +424 99 2 4 100 +425 76 29 91 30 +426 32 56 91 31 +427 55 37 107 38 +428 75 35 107 36 +429 4 99 100 84 +430 99 64 2 100 +431 90 34 47 33 +432 27 67 28 90 +433 88 108 62 82 +434 74 14 36 106 +435 87 93 81 61 +436 65 18 1 97 +437 3 97 15 85 +438 11 44 49 102 +439 12 69 45 102 +440 13 54 37 106 +441 14 13 10 106 +442 13 9 10 106 +443 84 100 4 46 +444 4 99 84 16 +445 17 64 2 99 +446 100 64 2 43 +447 34 5 87 61 +448 27 6 81 87 +449 92 1 3 97 +450 7 8 108 88 +451 38 62 7 108 +452 35 108 8 82 +453 22 93 5 61 +454 23 81 6 93 +455 3 92 85 26 +456 1 19 65 92 +457 92 65 1 97 +458 97 3 92 85 +459 5 93 87 61 +460 93 81 6 87 +461 108 62 7 88 +462 88 8 108 82 +463 86 103 12 41 +464 40 11 66 103 +465 14 74 42 105 +466 39 54 13 105 +467 6 5 93 87 +468 9 13 37 106 +469 10 36 14 106 +470 8 88 30 82 +471 31 7 62 88 +$EndElements diff --git a/meshes/msh/unitTetra.msh b/meshes/msh/unitTetra.msh index 9d9183d..1cb058e 100644 --- a/meshes/msh/unitTetra.msh +++ b/meshes/msh/unitTetra.msh @@ -50,7 +50,7 @@ $Nodes 3 1 0 0 $EndNodes $Elements -8 8 1 8 +9 9 1 9 0 1 15 1 1 1 0 2 15 1 @@ -67,4 +67,6 @@ $Elements 7 1 2 3 2 4 2 1 8 1 4 2 +3 1 4 1 +9 1 3 4 2 $EndElements diff --git a/poisson/CMakeLists.txt b/poisson/CMakeLists.txt index 0d2264f..1c01563 100644 --- a/poisson/CMakeLists.txt +++ b/poisson/CMakeLists.txt @@ -43,6 +43,7 @@ arcane_add_arcane_libraries_to_target(Poisson) target_include_directories(Poisson PUBLIC . ${CMAKE_CURRENT_BINARY_DIR}) configure_file(Poisson.config ${CMAKE_CURRENT_BINARY_DIR} COPYONLY) configure_file(Test.poisson.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONLY) +configure_file(Test.poisson.3D.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONLY) configure_file(Test.poisson.direct.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONLY) configure_file(Test.poisson.neumann.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONLY) configure_file(Test.poisson.trilinos.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONLY) @@ -51,6 +52,7 @@ configure_file(Test.poisson.hypre_direct.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONL configure_file(Test.poisson.petsc.arc ${CMAKE_CURRENT_BINARY_DIR} COPYONLY) configure_file(${MSH_DIR}/L-shape.msh ${CMAKE_CURRENT_BINARY_DIR} COPYONLY) configure_file(${MSH_DIR}/random.msh ${CMAKE_CURRENT_BINARY_DIR} COPYONLY) +configure_file(${MSH_DIR}/L-shape-3D.msh ${CMAKE_CURRENT_BINARY_DIR} COPYONLY) target_link_libraries(Poisson PUBLIC FemUtils) @@ -77,7 +79,7 @@ file(COPY "tests/" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) enable_testing() - +add_test(NAME [poisson]poisson_3D_Dirichlet COMMAND Poisson Test.poisson.3D.arc) add_test(NAME [poisson]poisson_direct COMMAND Poisson Test.poisson.direct.arc) if(FEMUTILS_HAS_SOLVER_BACKEND_TRILINOS) diff --git a/poisson/CsrBiliAssembly.hxx b/poisson/CsrBiliAssembly.hxx index 006e781..2c1bbe6 100644 --- a/poisson/CsrBiliAssembly.hxx +++ b/poisson/CsrBiliAssembly.hxx @@ -34,7 +34,12 @@ _buildMatrixCsr() } */ - Int32 nnz = nbFace() * 2 + nbNode(); + Int32 nnz; + if (options()->meshType == "TETRA4") + nnz = nbFace() * 18 + nbNode()*4; + else if (options()->meshType == "TRIA3") + nnz = nbFace() * 2 + nbNode(); + m_csr_matrix.initialize(m_dof_family, nnz, nbNode()); auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView()); //We iterate through the node, and we do not sort anymore : we assume the nodes ID are sorted, and we will iterate throught the column to avoid making < and > comparison @@ -159,4 +164,46 @@ _assembleCsrBilinearOperatorTRIA3() } /*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ \ No newline at end of file +/*---------------------------------------------------------------------------*/ + +void FemModule:: +_assembleCsrBilinearOperatorTETRA4() +{ + + + { + // Build the csr matrix + _buildMatrixCsr(); + } + + auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView()); + + ENUMERATE_ (Cell, icell, allCells()) { + Cell cell = *icell; + + FixedMatrix<4, 4> K_e; + { + K_e = _computeElementMatrixTETRA4(cell); // element stifness matrix + } + + //Timer::Action timer_action(m_time_stats, "CsrAddToGlobalMatrix"); + 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_csr_matrix.matrixAddValue(node_dof.dofId(node1, 0), node_dof.dofId(node2, 0), v); + } + ++n2_index; + } + ++n1_index; + } + + } +} + +/*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*/ diff --git a/poisson/FemModule.cc b/poisson/FemModule.cc index 4963647..ce88e21 100644 --- a/poisson/FemModule.cc +++ b/poisson/FemModule.cc @@ -231,7 +231,7 @@ compute() string_list.add("draw"); string_list.add("-draw_pause"); string_list.add("-10"); -*/ + */ CommandLineArguments args(string_list); m_linear_system.setSolverCommandLineArguments(args); } @@ -372,10 +372,44 @@ _doStationarySolve() _updateBoundayConditions(); // Assemble the FEM bilinear operator (LHS - matrix A) - if (options()->meshType == "QUAD4") + if (options()->meshType == "QUAD4"){ _assembleBilinearOperatorQUAD4(); + } + +/* + if (m_use_csr){ + m_linear_system.clearValues(); + _assembleCsrBilinearOperatorTETRA4(); + if (m_cache_warming != 1) { + m_time_stats->resetStats("AssembleCsrBilinearOperatorTria3"); + for (cache_index = 1; cache_index < m_cache_warming; cache_index++){ + m_linear_system.clearValues(); + _assembleCsrBilinearOperatorTETRA4(); + } + } + m_csr_matrix.translateToLinearSystem(m_linear_system); + } +*/ else { + if (m_use_legacy) { + m_linear_system.clearValues(); + if (options()->meshType == "TETRA4") + _assembleBilinearOperatorTETRA4(); + else if (options()->meshType == "TRIA3") + _assembleBilinearOperatorTRIA3(); + if (m_cache_warming != 1) { + m_time_stats->resetStats("AssembleLegacyBilinearOperatorTria3"); + for (cache_index = 1; cache_index < m_cache_warming; cache_index++) { + m_linear_system.clearValues(); + if (options()->meshType == "TETRA4") + _assembleBilinearOperatorTETRA4(); + else if (options()->meshType == "TRIA3") + _assembleBilinearOperatorTRIA3(); + } + } + } + #ifdef USE_CUSPARSE_ADD if (m_use_cusparse_add) { cusparseHandle_t handle; @@ -432,17 +466,6 @@ _doStationarySolve() } m_csr_matrix.translateToLinearSystem(m_linear_system); } - if (m_use_legacy) { - m_linear_system.clearValues(); - _assembleBilinearOperatorTRIA3(); - if (m_cache_warming != 1) { - m_time_stats->resetStats("AssembleLegacyBilinearOperatorTria3"); - for (cache_index = 1; cache_index < m_cache_warming; cache_index++) { - m_linear_system.clearValues(); - _assembleBilinearOperatorTRIA3(); - } - } - } #ifdef ARCANE_HAS_ACCELERATOR if (m_use_csr_gpu) { @@ -492,9 +515,9 @@ _doStationarySolve() m_csr_matrix.translateToLinearSystem(m_linear_system); _translateRhs(); } - else + else{ _assembleLinearOperator(); - + } _solve(); @@ -523,7 +546,7 @@ _getMaterialParameters() f = options()->f(); ElementNodes = 3.; - if (options()->meshType == "QUAD4") + if (options()->meshType == "QUAD4" || options()->meshType == "TETRA4") ElementNodes = 4.; } @@ -643,6 +666,9 @@ _checkCellType() if (options()->meshType == "QUAD4") { type = IT_Quad4; } + else if (options()->meshType == "TETRA4") { + type = IT_Tetraedron4; + } else { type = IT_Triangle3; } @@ -670,7 +696,7 @@ _updateBoundayConditions() void FemModule:: _assembleLinearOperator() { - info() << "Assembly of FEM linear operator "; + info() << "Assembly of FEM linear operator "; info() << "Applying Dirichlet boundary condition via penalty method "; // time registration @@ -843,14 +869,26 @@ _assembleLinearOperator() // $int_{Omega}(f*v^h)$ // only for noded that are non-Dirichlet //---------------------------------------------- - ENUMERATE_ (Cell, icell, allCells()) { - Cell cell = *icell; + if (options()->meshType == "TRIA3"){ + ENUMERATE_ (Cell, icell, allCells()) { + Cell cell = *icell; + Real area = _computeAreaTriangle3(cell); + for (Node node : cell.nodes()) { + if (!(m_u_dirichlet[node]) && node.isOwn()) { + rhs_values[node_dof.dofId(node, 0)] += f * area / ElementNodes; + } + } + } + } - Real area = _computeAreaTriangle3(cell); - for (Node node : cell.nodes()) { - if (!(m_u_dirichlet[node]) && node.isOwn()) { - // Original code - rhs_values[node_dof.dofId(node, 0)] += f * area / ElementNodes; + if (options()->meshType == "TETRA4"){ + ENUMERATE_ (Cell, icell, allCells()) { + Cell cell = *icell; + Real area = _computeAreaTetra4(cell); + for (Node node : cell.nodes()) { + if (!(m_u_dirichlet[node]) && node.isOwn()) { + rhs_values[node_dof.dofId(node, 0)] += f * area / ElementNodes; + } } } } @@ -1743,6 +1781,27 @@ _computeEdgeLength2(Face face) /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ + +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; +} + +/*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*/ + ARCCORE_HOST_DEVICE Real2 FemModule:: _computeEdgeNormal2Gpu(FaceLocalId iface, IndexedFaceNodeConnectivityView fnc, @@ -1879,6 +1938,111 @@ _computeElementMatrixQUAD4(Cell cell) /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ +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; +} + +/*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*/ + +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; + } + } + +} + +/*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*/ + void FemModule:: _assembleBilinearOperatorQUAD4() { diff --git a/poisson/FemModule.h b/poisson/FemModule.h index 72b7af7..a6ab324 100644 --- a/poisson/FemModule.h +++ b/poisson/FemModule.h @@ -248,6 +248,7 @@ class FemModule void _checkCellType(); void _assembleBilinearOperatorTRIA3(); void _assembleBilinearOperatorQUAD4(); + void _assembleBilinearOperatorTETRA4(); void _solve(); void _initBoundaryconditions(); void _assembleLinearOperator(); @@ -259,12 +260,14 @@ class FemModule void _benchBuildRow(); Real _readTimeFromJson(String main_time, String sub_time); 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); - //#ifdef ARCANE_HAS_ACCELERATOR + public: void _applyDirichletBoundaryConditionsGpu(); @@ -286,7 +289,7 @@ class FemModule private: - //#endif + #ifdef USE_CUSPARSE_ADD void printCsrMatrix(std::string fileName, cusparseCsr csr, bool is_coo); void _computeCusparseElementMatrix(cusparseCsr& result, cusparseCsr& global, Cell icell, cusparseHandle_t handle, IndexedNodeDoFConnectivityView node_dof, @@ -326,6 +329,7 @@ class FemModule #endif void _assembleCsrBilinearOperatorTRIA3(); + void _assembleCsrBilinearOperatorTETRA4(); void _buildMatrixCsr(); public: diff --git a/poisson/Test.poisson.3D.arc b/poisson/Test.poisson.3D.arc new file mode 100644 index 0000000..aba5674 --- /dev/null +++ b/poisson/Test.poisson.3D.arc @@ -0,0 +1,38 @@ + + + + Sample + PoissonLoop + + + + 1 + + U + + + + + + L-shape-3D.msh + + + + + true + 1.0 + test_3D_L-shape_poisson.txt + TETRA4 + + bot + 50.0 + + + bc + 10.0 + + + petsc + + + diff --git a/poisson/tests/test_3D_L-shape_poisson.txt b/poisson/tests/test_3D_L-shape_poisson.txt new file mode 100644 index 0000000..bf6397b --- /dev/null +++ b/poisson/tests/test_3D_L-shape_poisson.txt @@ -0,0 +1,108 @@ +1 50 +2 50 +3 50 +4 50 +5 108089.067464837 +6 108089.368096104 +7 84544.0611838288 +8 84612.1578854208 +9 10 +10 10 +11 81932.2855675261 +12 82528.0563642635 +13 10 +14 10 +15 50 +16 50 +17 50 +18 50 +19 40195.4718795868 +20 71152.3155106399 +21 91431.8303087035 +22 103711.128094321 +23 103651.193098788 +24 91693.0490337547 +25 70772.0478435848 +26 40556.6361500248 +27 106794.570914974 +28 103194.404817591 +29 96753.7560372044 +30 89511.8979601444 +31 88635.9714930756 +32 96693.9294356658 +33 103411.590661594 +34 106984.310218392 +35 77933.251009223 +36 54226.5784548942 +37 54166.6216888929 +38 77646.316957521 +39 57482.3692176502 +40 78246.1316590261 +41 78110.2458079007 +42 57349.4249578766 +43 34246.6950269383 +44 61900.1568450923 +45 61842.3185014287 +46 34243.9197095375 +47 103349.138611966 +48 74149.5201401205 +49 71906.1628916932 +50 30182.1712007777 +51 88406.584464632 +52 91498.0522071363 +53 82236.7491242031 +54 53041.7543454436 +55 72928.7223373876 +56 89604.7408391188 +57 55338.0206474763 +58 98316.1518446525 +59 98770.2667805917 +60 48052.0630627939 +61 105461.001447101 +62 83304.7709487885 +63 97408.1352050961 +64 25855.5406563433 +65 26950.485365453 +66 84386.8320677527 +67 103152.728804871 +68 73465.0679874443 +69 71962.85079717 +70 30150.2467830329 +71 88710.5298321351 +72 91526.7667330718 +73 81756.9150910979 +74 53753.778622362 +75 73028.8861227429 +76 88894.6169693091 +77 55425.6224295762 +78 98524.3037809767 +79 98847.110848987 +80 47854.0677828412 +81 105374.990794904 +82 84022.7208603989 +83 97281.5819081245 +84 25835.8016660597 +85 27000.9065713509 +86 84014.1889111783 +87 107681.303635058 +88 86061.3293681357 +89 100879.803773252 +90 105847.463899652 +91 93089.2672194831 +92 21371.1822745454 +93 106682.178278461 +94 82819.6504668915 +95 57583.0381058618 +96 98766.536582697 +97 50 +98 50 +99 50 +100 21268.1781940651 +101 48838.7403238286 +102 71359.4349424098 +103 80022.7446737088 +104 69281.6954098975 +105 43988.4162283485 +106 35316.7707752914 +107 68326.6395839078 +108 81693.9670839977