diff --git a/srcs/BEM/functionsBEM.cpp b/srcs/BEM/functionsBEM.cpp index 38a77d269d244ce8760cbc21438c68fe6c190f28..14ed1b34103e1d9398ffec64c9f39871db9242fc 100644 --- a/srcs/BEM/functionsBEM.cpp +++ b/srcs/BEM/functionsBEM.cpp @@ -32,12 +32,11 @@ void dispValuesFun(const std::vector<elementStruct> &elementVector) } } - -void getBC(const std::vector<double> &allCoord, std::vector<elementStruct> &elementVector, int &dirichletCount, const bool dispHierarchy) +void getBC(std::map<std::size_t, int> allNodeIndex, const std::vector<double> &allCoord, std::vector<elementStruct> &elementVector, int &dirichletCount, const bool dispHierarchy) { // Gets the dim and the tags of the physical groups gmsh::vectorpair physicalGroupsDimTags; - gmsh::model::getPhysicalGroups(physicalGroupsDimTags); + gmsh::model::getPhysicalGroups(physicalGroupsDimTags, 1); // Gets the boundary conditions from onelab and create a map: "name of the physical group" -> ("dirichlet"/"neumann", value) std::map<std::string, std::pair<std::string, double>> bcMap; @@ -67,6 +66,9 @@ void getBC(const std::vector<double> &allCoord, std::vector<elementStruct> &elem std::string name; gmsh::model::getPhysicalName(physicalGroupsDim, physicalGroupsTag, name); + if(bcMap.find(name) == bcMap.end()) + continue; + // Retrieves the type and the value of the boundary condition of this physical group. // /!\ Should check that this physical groups was well defined in the list of "Boundary Condition" ! std::string bcType = bcMap[name].first; @@ -85,26 +87,27 @@ void getBC(const std::vector<double> &allCoord, std::vector<elementStruct> &elem if(dispHierarchy) std::cout << "\t contains entity '" << entitiesTags[k] << "' "; - std::vector<std::size_t> elementTags; - std::vector<std::size_t> nodeTags; - gmsh::model::mesh::getElementsByType(1, elementTags, nodeTags, entitiesTags[k]); - + + std::vector<std::size_t> elementTags1; + std::vector<std::size_t> nodeTags1; + std::vector<std::size_t> elementTags8; + std::vector<std::size_t> nodeTags8; + gmsh::model::mesh::getElementsByType(1, elementTags1, nodeTags1, entitiesTags[k]); + gmsh::model::mesh::getElementsByType(8, elementTags8, nodeTags8, entitiesTags[k]); + // Print the number of element of the given entity (if dispHierarchy != 0). if(dispHierarchy) - std::cout << "which has " << elementTags.size() << " elements:\n"; + std::cout << "which has " << elementTags1.size() + elementTags8.size() << " elements:\n"; - for(int el = 0; el < elementTags.size(); el++) + for(int el = 0; el < elementTags1.size(); el++) { - // Stores the tag, the boundary type, the boundary value and the node tags of the element. - int elTag = elementTags[el]; + int elTag = elementTags1[el]; elementStruct element; element.tag = elTag; - element.bcType = bcType; - element.nodeTags[0] = nodeTags[2*el]; - element.nodeTags[1] = nodeTags[2*el + 1]; - computeGeometricParam(element, allCoord); - // The elementVector should be "sorted" => push from front or from back depending on element type. + element.nodeTags[0] = nodeTags1[2*el]; + element.nodeTags[1] = nodeTags1[2*el + 1]; + computeGeometricParam(element, allNodeIndex, allCoord); if(bcType == "dirichlet") { element.bcValue[0] = bcValue; @@ -116,13 +119,41 @@ void getBC(const std::vector<double> &allCoord, std::vector<elementStruct> &elem element.bcValue[1] = bcValue; elementVector.push_back(element); } + // Print the tag of a given element and the tags of the nodes that composes it as well as their coordinates (if dispHierarchy != 0). + if(dispHierarchy) + { + std::cout << "\t\t-Element '" << elTag << "' (" << element.bcType << "):\n"; + std::cout << "\t\t\t-> Node '" << element.nodeTags[0] << "' (" << allCoord[3*allNodeIndex[element.nodeTags[0]]] << "; " << allCoord[3*allNodeIndex[element.nodeTags[0]] + 1] << ")\n"; + std::cout << "\t\t\t-> Node '" << element.nodeTags[1] << "' (" << allCoord[3*allNodeIndex[element.nodeTags[1]]] << "; " << allCoord[3*allNodeIndex[element.nodeTags[1]] + 1] << ")\n"; + } + } + for(int el = 0; el < elementTags8.size(); el++) + { + int elTag = elementTags8[el]; + elementStruct element; + element.tag = elTag; + element.nodeTags[0] = nodeTags8[3*el]; + element.nodeTags[1] = nodeTags8[3*el + 1]; + computeGeometricParam(element, allNodeIndex, allCoord); + if(bcType == "dirichlet") + { + element.bcValue[0] = bcValue; + elementVector.insert(elementVector.begin(), element); + dirichletCount++; + } + else if(bcType == "neumann") + { + element.bcValue[1] = bcValue; + elementVector.push_back(element); + } + // Print the tag of a given element and the tags of the nodes that composes it as well as their coordinates (if dispHierarchy != 0). if(dispHierarchy) { std::cout << "\t\t-Element '" << elTag << "' (" << element.bcType << "):\n"; - std::cout << "\t\t\t-> Node '" << element.nodeTags[0] << "' (" << allCoord[3*(element.nodeTags[0]-1)] << "; " << allCoord[3*(element.nodeTags[0]-1) + 1] << ")\n"; - std::cout << "\t\t\t-> Node '" << element.nodeTags[1] << "' (" << allCoord[3*(element.nodeTags[1]-1)] << "; " << allCoord[3*(element.nodeTags[1]-1) + 1] << ")\n"; + std::cout << "\t\t\t-> Node '" << element.nodeTags[0] << "' (" << allCoord[3*allNodeIndex[element.nodeTags[0]]] << "; " << allCoord[3*allNodeIndex[element.nodeTags[0]] + 1] << ")\n"; + std::cout << "\t\t\t-> Node '" << element.nodeTags[1] << "' (" << allCoord[3*allNodeIndex[element.nodeTags[1]]] << "; " << allCoord[3*allNodeIndex[element.nodeTags[1]] + 1] << ")\n"; } } if(dispHierarchy) @@ -132,16 +163,19 @@ void getBC(const std::vector<double> &allCoord, std::vector<elementStruct> &elem } -void computeGeometricParam(elementStruct &element, const std::vector<double> &allCoord) +void computeGeometricParam(elementStruct &element, std::map<std::size_t, int> allNodeIndex, const std::vector<double> &allCoord) { int i1Tag = element.nodeTags[0]; int i2Tag = element.nodeTags[1]; + int i1Index = allNodeIndex[i1Tag]; + int i2Index = allNodeIndex[i2Tag]; + // Get the coordinates of the nodes from their tag. - element.x1 = allCoord[3*(i1Tag-1)]; - element.y1 = allCoord[3*(i1Tag-1) + 1]; - element.x2 = allCoord[3*(i2Tag-1)]; - element.y2 = allCoord[3*(i2Tag-1) + 1]; + element.x1 = allCoord[3*i1Index]; + element.y1 = allCoord[3*i1Index + 1]; + element.x2 = allCoord[3*i2Index]; + element.y2 = allCoord[3*i2Index + 1]; element.midX = (element.x1 + element.x2) / 2; element.midY = (element.y1 + element.y2) / 2; diff --git a/srcs/BEM/functionsBEM.hpp b/srcs/BEM/functionsBEM.hpp index 5b0f375ded4ec865e34b3b9337abe0c800043a5c..bc9276850678d916f6a5a6dee5c4e643a765474b 100644 --- a/srcs/BEM/functionsBEM.hpp +++ b/srcs/BEM/functionsBEM.hpp @@ -1,5 +1,6 @@ #include <iostream> #include <vector> // Why should I include this ? +#include <map> #include <Eigen/Dense> // Each element has all its information gathered as in this structure. @@ -7,7 +8,7 @@ struct elementStruct{ int tag; // tag of the element std::string bcType; // "dirichlet" or "neumann" (at the end, should replace this by an int or boolean) double bcValue [2]; // values of [u, u_n] - int nodeTags [2]; // tags of the [first, second] nodes of the element + std::size_t nodeTags [2]; // tags of the [first, second] nodes of the element double length; double x1, y1, x2, y2; @@ -19,8 +20,8 @@ void dispNodeCoordFun(const std::vector<size_t> &allNodeTags, const std::vector< void dispMatricesFun(const Eigen::MatrixXd &A, const Eigen::MatrixXd &B, const Eigen::MatrixXd &c, const Eigen::MatrixXd &b, const Eigen::MatrixXd &x); void dispValuesFun(const std::vector<elementStruct> &elementVector); -void getBC(const std::vector<double> &allCoord, std::vector<elementStruct> &elementVector, int &dirichletCount, const bool dispHierarchy); -void computeGeometricParam(elementStruct &element, const std::vector<double> &allCoord); +void getBC(std::map<std::size_t, int> allNodeIndex, const std::vector<double> &allCoord, std::vector<elementStruct> &elementVector, int &dirichletCount, const bool dispHierarchy); +void computeGeometricParam(elementStruct &element, std::map<std::size_t, int> allNodeIndex, const std::vector<double> &allCoord); double computeH(const double x, const double y, const elementStruct &element); double computeG(const double x, const double y, const elementStruct &element, const std::string integrationType); diff --git a/srcs/BEM/solverBEM.cpp b/srcs/BEM/solverBEM.cpp index f0d17049b0fa94534cbc5508d791986d9a3fadff..02c8bfdfd6b9bcb289442228436734453af09cc0 100644 --- a/srcs/BEM/solverBEM.cpp +++ b/srcs/BEM/solverBEM.cpp @@ -25,17 +25,20 @@ int solverBEM(int argc, char **argv, std::vector<size_t> &finalElementTags, std: std::vector<elementStruct> elementVector; // Gets the tags, the coordinates and the parametric coordinates of the nodes. + std::map<std::size_t, int> allNodeIndex; // nodeTag -> position in allNodeTags and allCoord std::vector<std::size_t> allNodeTags; std::vector<double> allCoord; std::vector<double> allParametricCoord; gmsh::model::mesh::getNodes(allNodeTags, allCoord, allParametricCoord); + for(int i = 0; i < allNodeTags.size(); i++) + allNodeIndex[allNodeTags[i]] = i; // Print the coordinates of all nodes if dispNodeCoord != 0 if(dispNodeCoord) dispNodeCoordFun(allNodeTags, allCoord); int dirichletCount = 0; - getBC(allCoord, elementVector, dirichletCount, dispHierarchy); + getBC(allNodeIndex, allCoord, elementVector, dirichletCount, dispHierarchy); // Declaration of matrices A, B and vector c (the matrix type is perhaps not the most efficient). Eigen::MatrixXd A = Eigen::MatrixXd::Zero(elementVector.size(), elementVector.size()); diff --git a/srcs/complex_validation.geo b/srcs/complex_validation.geo new file mode 100644 index 0000000000000000000000000000000000000000..6d7cb994b6a3044f89a63cbe842bc619f96527d1 --- /dev/null +++ b/srcs/complex_validation.geo @@ -0,0 +1,51 @@ +//THIS CODE DEFINES THE GEOMETRY OF A CLAMPED BAR WITH TWO TYPES OF MESHES, TRIANGULAR AND QUADRANGULAR + +Lx = 4; +Ly = 2; +nx = 40; +ny = 20; + +Point(1) = {0, 0, 0, 0.1}; +Point(2) = {Lx, 0, 0, 0.2}; +Point(3) = {Lx, Ly, 0, 0.2}; +Point(4) = {0, Ly, 0, 0.1}; +Point(5) = {Lx/2, 0, 0, 0.15}; +Point(6) = {Lx/2, Ly, 0, 0.15}; + +Line(1) = {1, 5}; +Line(2) = {5, 2}; +Line(3) = {2, 3}; +Line(4) = {3, 6}; +Line(5) = {6, 4}; +Line(6) = {4, 1}; +Line(7) = {5, 6}; +Line(8) = {6, 5}; + +Curve Loop(1) = {1, 7, 5, 6}; +Curve Loop(2) = {2, 3, 4, 8}; + +Plane Surface(1) = {1}; +Plane Surface(2) = {2}; + +Transfinite Curve {1, 5} = nx+1 Using Progression 1; +Transfinite Curve {7, 6} = ny+1 Using Progression 1; +Transfinite Curve {2, 4} = nx+1 Using Progression 1; +Transfinite Curve {3, 8} = ny+1 Using Progression 1; +Transfinite Surface {1}; +Transfinite Surface {2}; + + +Mesh.ElementOrder = 2; + +Recombine Surface {1}; // quads instead of triangles + +Physical Curve("left_edge", 5) = {6}; +Physical Curve("right_edge", 6) = {3}; +Physical Curve("bottom_edge", 7) = {1, 2}; +Physical Curve("top_edge", 8) = {4, 5}; +Physical Curve("between", 9) = {7, 8}; + +SetNumber("Boundary Conditions/left_edge/dirichlet", 200.); +SetNumber("Boundary Conditions/right_edge/dirichlet", 100.); +SetNumber("Boundary Conditions/bottom_edge/neumann", 0.); +SetNumber("Boundary Conditions/top_edge/neumann", 0.); \ No newline at end of file