From 4db1b11343382c387a4183235caf2ebc1dfd67d4 Mon Sep 17 00:00:00 2001
From: Pierre-Yves <PYGoffin@student.uliege.be>
Date: Wed, 30 Mar 2022 23:35:55 +0200
Subject: [PATCH] Update of the BEM code in order to be able to deal with
multiple meshes and order 2 elements. The code is however not able to deal
with a line inside the domain that has been assigned a boundary condition.
---
srcs/BEM/functionsBEM.cpp | 80 ++++++++++++++++++++++++++-----------
srcs/BEM/functionsBEM.hpp | 7 ++--
srcs/BEM/solverBEM.cpp | 5 ++-
srcs/complex_validation.geo | 51 +++++++++++++++++++++++
4 files changed, 116 insertions(+), 27 deletions(-)
create mode 100644 srcs/complex_validation.geo
diff --git a/srcs/BEM/functionsBEM.cpp b/srcs/BEM/functionsBEM.cpp
index 38a77d2..14ed1b3 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 5b0f375..bc92768 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 f0d1704..02c8bfd 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 0000000..6d7cb99
--- /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
--
GitLab