Skip to content
Snippets Groups Projects
Commit 4db1b113 authored by Goffin Pierre-Yves's avatar Goffin Pierre-Yves
Browse files

Update of the BEM code in order to be able to deal with multiple meshes and... 

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.
parent 13e61391
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
srcs/BEM/functionsBEM.cpp
+ 57
− 23
View file @ 4db1b113
......@@ -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;
......
srcs/BEM/functionsBEM.hpp
+ 4
− 3
View file @ 4db1b113
#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);
......
srcs/BEM/solverBEM.cpp
+ 4
− 1
View file @ 4db1b113
......@@ -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());
......
srcs/complex_validation.geo 0 → 100644
+ 51
− 0
View file @ 4db1b113
//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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment