Skip to content
Snippets Groups Projects

Py bem

Closed Goffin Pierre-Yves requested to merge PY_BEM into master
11 files
+ 2037
1558
Compare changes
  • Side-by-side
  • Inline
Files
11
+ 311
0
#include <gmsh.h>
#include <iostream>
#include <map>
#include <sstream>
#include <Eigen/Dense>
#include <cmath>
#include <iomanip>
#include <stdio.h>
#include <chrono>
#include "functionsBEM.hpp"
#include <omp.h>
#include <algorithm>
#include <list>
// Fills the "boundaryConditions" map associating to each physical group (a string) the type of boundary condition ('true':
// dirichlet, 'false': neumann) and its value (if any). "domainName" is the name of the current BEM domain.
void getBC(std::map<std::string, std::pair<bool, double>> &boundaryConditions,
const std::string domainName)
{
// Extracts all the boundary conditions from "ONELAB".
std::vector<std::string> bcKeys;
gmsh::onelab::getNames(bcKeys, "Boundary Conditions.+");
for(auto &key : bcKeys)
{
std::vector<double> values;
gmsh::onelab::getNumber(key, values);
std::stringstream ss(key);
std::vector<std::string> words;
std::string word;
while(getline(ss, word, '/'))
words.push_back(word);
// Only focus on the right BEM domain.
if(words.size() == 4 && words[2].compare(domainName) == 0)
{
if(words[3].compare("dirichlet") == 0)
boundaryConditions[words[1]] = {true, values[0]};
else if(words[3].compare("neumann") == 0)
boundaryConditions[words[1]] = {false, values[0]};
}
}
}
// Fills the "elementVector" vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of the
// current BEM domain, with the appropriate information for each element (the bcValue will be computed later). Also fills the
// "physicalGroups" map associating to each (physical group) name, its tag and dimension. "domainName" is the name of the current
// BEM domain. "nodalDisplacements" is a map associating the (x,y) displacement to a serie of node tags. Returns the number of
// 'dirichlet elements' (useful for the construction/resolution of the system).
std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
std::map<std::string, std::pair<int, int>> &physicalGroups,
std::map<int, std::pair<double, double>> &nodalDisplacements,
const std::string domainName)
{
std::size_t nbDirichletElements = 0; // Will be incremented later
std::map<std::string, std::pair<bool, double>> boundaryConditions;
getBC(boundaryConditions, domainName);
// Gets the dim and the tags of the physical groups.
gmsh::vectorpair physicalGroupsDimTags;
gmsh::model::getPhysicalGroups(physicalGroupsDimTags);
for(std::size_t i = 0; i < physicalGroupsDimTags.size(); i++)
{
// Fills the "physicalGroups" map that associates the name to the {dim, tag}.
int physicalGroupDim = physicalGroupsDimTags[i].first;
int physicalGroupTag = physicalGroupsDimTags[i].second;
std::string physicalGroupName;
gmsh::model::getPhysicalName(physicalGroupDim, physicalGroupTag, physicalGroupName);
physicalGroups[physicalGroupName] = {physicalGroupDim, physicalGroupTag};
gmsh::model::getPhysicalName(physicalGroupDim, physicalGroupTag, physicalGroupName);
// If no boundary condition on a given physical curve, don't do anything.
if(boundaryConditions.find(physicalGroupName) == boundaryConditions.end())
continue;
// Retrieves the type and the value of the boundary condition of this physical group.
bool bcType = boundaryConditions[physicalGroupName].first;
double bcValue = boundaryConditions[physicalGroupName].second;
// Loops over the tags of the entities and the types of the elements for retrieving the relevant information (element
// tags and nodes tags) in gmsh.
std::vector<int> entityTags;
gmsh::model::getEntitiesForPhysicalGroup(physicalGroupDim, physicalGroupTag, entityTags);
for (std::size_t i = 0; i < entityTags.size(); i++)
{
std::vector<int> elementTypes = {1, 8}; // Hardcoded types of 2/3-nodes 1D elements.
for(std::size_t j = 0; j < elementTypes.size(); j++)
{
int nbNodesPerElement = 2;
if(elementTypes[j] == 8)
nbNodesPerElement = 3;
// Get the element tags and the node tags of each element.
std::vector<std::size_t> elementTags;
std::vector<std::size_t> elementNodeTags;
gmsh::model::mesh::getElementsByType(elementTypes[j], elementTags, elementNodeTags, entityTags[i]);
// Get the coordinates of each node (extracted above).
std::vector<std::size_t> nodeTags;
std::vector<double> nodeCoords;
std::vector<double> parametricNodeCoords;
gmsh::model::mesh::getNodesByElementType(elementTypes[j], nodeTags, nodeCoords, parametricNodeCoords, entityTags[i], false);
// Constructs a map associating the index of each node tag in the "nodeCoords" vector.
std::map<std::size_t, std::size_t> nodeIndices;
for(std::size_t k = 0; k < nodeTags.size(); k++)
nodeIndices[nodeTags[k]] = k;
// Computes each element data and fills the "elementVector" vector.
#pragma omp parallel shared(nbDirichletElements, elementTags, nodeTags, nodalDisplacements, nodeIndices, nodeCoords, bcType, bcValue)
{
std::vector<elementStruct> localDirichletElementVector;
std::vector<elementStruct> localNeumannElementVector;
// Computes the element data
#pragma omp for schedule(guided)
for(std::size_t el = 0; el < elementTags.size(); el++)
{
std::size_t elTag = elementTags[el];
elementStruct element;
element.tag = elTag;
element.nodeTags[0] = nodeTags[nbNodesPerElement*el];
element.nodeTags[1] = nodeTags[nbNodesPerElement*el + 1];
if(nbNodesPerElement == 3)
element.midNodeTag = nodeTags[nbNodesPerElement*el + 2];
else
element.midNodeTag = 0; // "midNodeTag" set to '0' (default value) if no mid node on the element.
std::pair<double, double> nodalDisp1 = {0.0, 0.0}, nodalDisp2 = {0.0, 0.0};
if(nodalDisplacements.count(element.nodeTags[0]))
nodalDisp1 = nodalDisplacements[element.nodeTags[0]];
if(nodalDisplacements.count(element.nodeTags[1]))
nodalDisp2 = nodalDisplacements[element.nodeTags[1]];
computeGeometricParam(element, nodeIndices, nodeCoords, nodalDisp1, nodalDisp2);
// Pushes the element at the back of lists (private for each thread)/
if(bcType)
{
element.bcValue[0] = bcValue;
localDirichletElementVector.push_back(element);
}
else
{
element.bcValue[1] = bcValue;
localNeumannElementVector.push_back(element);
}
}
// Add the elements to the "elementVector" vector. It is important that the 'dirichlet elements' are first in
// the vector.
#pragma omp critical
{
elementVector.insert(elementVector.begin(), localDirichletElementVector.begin(), localDirichletElementVector.end());
elementVector.insert(elementVector.end(), localNeumannElementVector.begin(), localNeumannElementVector.end());
nbDirichletElements += localDirichletElementVector.size();
}
}
}
}
}
return nbDirichletElements;
}
// Computes the useful data of a given element and stores it in "element". "nodeIndinces" is a map associating the index of any
// node tag for accessing its coordinates in the "nodeCoords" vector. "nodalDisp1" and "nodalDisp2" are pairs representing the
// potential displacement of the first and second node (respectivelly) of the element. If an element has 3 nodes, the first and
// second nodes are considered to be the extremities of that element.
void computeGeometricParam(elementStruct &element,
std::map<std::size_t, std::size_t> &nodeIndices,
const std::vector<double> &nodeCoords,
const std::pair<double, double> &nodalDisp1,
const std::pair<double, double> &nodalDisp2)
{
std::size_t i1Tag = element.nodeTags[0];
std::size_t i2Tag = element.nodeTags[1];
std::size_t i1Index = nodeIndices[i1Tag];
std::size_t i2Index = nodeIndices[i2Tag];
// Get the coordinates of the nodes from their tag.
element.x1 = nodeCoords[3*i1Index] + nodalDisp1.first;
element.x2 = nodeCoords[3*i2Index] + nodalDisp2.first;
element.y1 = nodeCoords[3*i1Index + 1] + nodalDisp1.second;
element.y2 = nodeCoords[3*i2Index + 1] + nodalDisp2.second;
element.midX = (element.x1 + element.x2) / 2;
element.midY = (element.y1 + element.y2) / 2;
element.deltaX = element.x2 - element.x1;
element.deltaY = element.y2 - element.y1;
element.length = sqrt(element.deltaX * element.deltaX + element.deltaY * element.deltaY);
}
// Fills the "A", "B", "c" matrices (/vector) that represent the linear system. "elementVector" is a vector storing the boundary
// elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM domain. "nbDirichletElements" is the number of
// 'dirichlet elements'. "solutionType" indicates whether the elements of the matrices are computed analytically ('true') or
// numerically ('false'). "integrationType" is the type of numerical integration used (in case of numerical computation).
void fillSystem(Eigen::MatrixXd &A,
Eigen::MatrixXd &B,
Eigen::MatrixXd &c,
const std::vector<elementStruct> &elementVector,
const std::size_t nbDirichletElements,
const std::string integrationType,
const bool solutionType)
{
double h;
double g;
// Double loop on the elements for computing each elements of the linear system.
#pragma omp parallel for collapse(2) schedule(guided) shared(elementVector, A, B, c) private(h, g)
for(std::size_t i = 0; i < elementVector.size(); i++)
{
for(std::size_t j = 0; j < elementVector.size(); j++)
{
// Will be computed and stored in A and B.
h = 0;
g = 0;
if(i == j) // Formula for diagonal elements.
{
h = -0.5;
g = (elementVector[i].length/(2*M_PI)) * (log(elementVector[i].length/2)-1);
}
else // General case
{
h = computeAnalyticalH(elementVector[i].midX, elementVector[i].midY, elementVector[j]);
if(!solutionType)
g = computeNumericalG(elementVector[i].midX, elementVector[i].midY, elementVector[j], integrationType);
else
g = computeAnalyticalG(elementVector[i].midX, elementVector[i].midY, elementVector[j]);
}
// Storage of h and g at the right place (in matrices A and B and vector c).
if(j < nbDirichletElements)
{
A(i,j) = g;
B(i,j) = h;
c(j) = elementVector[j].bcValue[0];
}
else
{
A(i,j) = -h;
B(i,j) = -g;
c(j) = elementVector[j].bcValue[1];
}
}
}
}
// Fills the "electrostaticPressure" map associating each (boundary) element tag with its electrostatic pressure. Also sets the
// unknown of each element (either the potential of the electric field) of the "elementVector" vector storing the boundary
// elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM domain. "x" is the solution of the linear
// system. "nbDirichletElements" is the number of 'dirichlet elements', "domainName" is the name of the considered BEM domain.
void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
std::vector<elementStruct> &elementVector,
const Eigen::MatrixXd &x,
const std::size_t nbDirichletElements,
const std::string domainName)
{
// Retrieves the dielectric permittivity of the BEM domain. If not set in the .geo file, set to the dielectric permittivity
// of vacuum (with a warning).
std::vector<std::string> keys;
gmsh::onelab::getNames(keys, "(Materials).+");
double epsilon = 8.8541878128e-12; // dielectric permittivity
bool setEpsilon = false;
for (auto &key : keys)
{
std::vector<double> value;
gmsh::onelab::getNumber(key, value);
// Expected key structure is "type/group_name/field"
std::stringstream ss(key);
std::vector<std::string> words;
std::string word;
while(std::getline(ss, word, '/')) // read string until '/'
words.push_back(word);
if(words.size()==3)
{
if(words[2].compare("Epsilon") == 0 && words[1].compare(domainName) == 0)
{
epsilon = value[0];
setEpsilon = true;
}
}
}
if(!setEpsilon)
std::cout << "\tWarning: You forgot to set the dielectric permittivity of the material for the domain '" << domainName << "' in the .geo file and it thus has been fixed to the dielectric permittivity of the vacuum: " << epsilon << " [F/m].\n";
#pragma omp parallel for schedule(guided)
for(size_t i = 0; i < elementVector.size(); i++)
{
// Stores the unknown value (potential or electric field) in the appropriate elements.
if(i < nbDirichletElements)
elementVector[i].bcValue[1] = x(i);
else
elementVector[i].bcValue[0] = x(i);
// Computes the electrostatic pressure and fills the map.
double pressure = epsilon * elementVector[i].bcValue[1] * elementVector[i].bcValue[1] / 2; // don't forget factor 2
#pragma omp critical
electrostaticPressure[elementVector[i].tag] = pressure;
}
}
\ No newline at end of file
Loading