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

comments BEM

parent 9736bb32
No related branches found
No related tags found
1 merge request!16Py bem
Expand all Hide whitespace changes
Inline Side-by-side
srcs/BEM/boundaryElements.cpp
+ 23
− 24
View file @ 81ad07f7
#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 "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': // 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. // dirichlet, 'false': neumann) and its value (if any). "domainName" is the name of the current BEM domain.
...@@ -44,12 +32,13 @@ void getBC(std::map<std::string, std::pair<bool, double>> &boundaryConditions, ...@@ -44,12 +32,13 @@ void getBC(std::map<std::string, std::pair<bool, double>> &boundaryConditions,
} }
// Fills the "elementVector" vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of the // 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 // current BEM domain, with the appropriate information for each element (the bcValue will be computed later). Also returns
// "physicalGroups" map associating to each (physical group) name, its tag and dimension. "domainName" is the name of the current // "domainPhysicalGroupDim" and "domainPhysicalGroupTag", the dimension and the tag (respectivelly) of the BEM domain of
// BEM domain. "nodalDisplacements" is a map associating the (x,y) displacement to a serie of node tags. Returns the number of // interest. "domainName" is the name of the current BEM domain. "nodalDisplacements" is a map associating the (x,y) displacement
// 'dirichlet elements' (useful for the construction/resolution of the system). // 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::size_t fillElementVector(std::vector<elementStruct> &elementVector,
std::map<std::string, std::pair<int, int>> &physicalGroups, int &domainPhysicalGroupDim,
int &domainPhysicalGroupTag,
std::map<int, std::pair<double, double>> &nodalDisplacements, std::map<int, std::pair<double, double>> &nodalDisplacements,
const std::string domainName) const std::string domainName)
{ {
...@@ -68,11 +57,14 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector, ...@@ -68,11 +57,14 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
int physicalGroupDim = physicalGroupsDimTags[i].first; int physicalGroupDim = physicalGroupsDimTags[i].first;
int physicalGroupTag = physicalGroupsDimTags[i].second; int physicalGroupTag = physicalGroupsDimTags[i].second;
std::string physicalGroupName; std::string physicalGroupName;
gmsh::model::getPhysicalName(physicalGroupDim, physicalGroupTag, physicalGroupName);
physicalGroups[physicalGroupName] = {physicalGroupDim, physicalGroupTag};
gmsh::model::getPhysicalName(physicalGroupDim, physicalGroupTag, physicalGroupName); gmsh::model::getPhysicalName(physicalGroupDim, physicalGroupTag, physicalGroupName);
if(physicalGroupName.compare(domainName) == 0)
{
domainPhysicalGroupDim = physicalGroupDim;
domainPhysicalGroupTag = physicalGroupTag;
}
// If no boundary condition on a given physical curve, don't do anything. // If no boundary condition on a given physical curve, don't do anything.
if(boundaryConditions.find(physicalGroupName) == boundaryConditions.end()) if(boundaryConditions.find(physicalGroupName) == boundaryConditions.end())
continue; continue;
...@@ -103,7 +95,8 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector, ...@@ -103,7 +95,8 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
std::vector<std::size_t> nodeTags; std::vector<std::size_t> nodeTags;
std::vector<double> nodeCoords; std::vector<double> nodeCoords;
std::vector<double> parametricNodeCoords; std::vector<double> parametricNodeCoords;
gmsh::model::mesh::getNodesByElementType(elementTypes[j], nodeTags, nodeCoords, parametricNodeCoords, entityTags[i], false); 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. // Constructs a map associating the index of each node tag in the "nodeCoords" vector.
std::map<std::size_t, std::size_t> nodeIndices; std::map<std::size_t, std::size_t> nodeIndices;
...@@ -156,8 +149,10 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector, ...@@ -156,8 +149,10 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
// the vector. // the vector.
#pragma omp critical #pragma omp critical
{ {
elementVector.insert(elementVector.begin(), localDirichletElementVector.begin(), localDirichletElementVector.end()); elementVector.insert(elementVector.begin(), localDirichletElementVector.begin(),
elementVector.insert(elementVector.end(), localNeumannElementVector.begin(), localNeumannElementVector.end()); localDirichletElementVector.end());
elementVector.insert(elementVector.end(), localNeumannElementVector.begin(),
localNeumannElementVector.end());
nbDirichletElements += localDirichletElementVector.size(); nbDirichletElements += localDirichletElementVector.size();
} }
} }
...@@ -291,7 +286,11 @@ void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure, ...@@ -291,7 +286,11 @@ void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
} }
} }
if(!setEpsilon) 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"; {
std::cout << "\tWarning: You forgot to set the dielectric permittivity of the material for the domain '" << domainName;
std::cout << "' in the .geo file and it thus has been fixed to the dielectric permittivity of the vacuum: " << epsilon;
std::cout << " [F/m].\n";
}
#pragma omp parallel for schedule(guided) #pragma omp parallel for schedule(guided)
for(size_t i = 0; i < elementVector.size(); i++) for(size_t i = 0; i < elementVector.size(); i++)
......
srcs/BEM/computations.cpp
+ 27
− 13
View file @ 81ad07f7
#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 "functionsBEM.hpp"
#include <omp.h>
#include <algorithm>
#include <list>
// Computes the 'h' element analytically. "x" and "y" are the coordinates at which 'h' should be computed. "element" is the
// element 'j' while will helps to compute 'h'.
double computeAnalyticalH(const double x, double computeAnalyticalH(const double x,
const double y, const double y,
const elementStruct &element) const elementStruct &element)
...@@ -41,6 +31,8 @@ double computeAnalyticalH(const double x, ...@@ -41,6 +31,8 @@ double computeAnalyticalH(const double x,
return alpha / (2*M_PI); return alpha / (2*M_PI);
} }
// Computes the 'g' element analytically. "x" and "y" are the coordinates at which 'g' should be computed. "element" is the
// element 'j' while will helps to compute 'g'.
double computeAnalyticalG(const double x, double computeAnalyticalG(const double x,
const double y, const double y,
const elementStruct &element) const elementStruct &element)
...@@ -58,6 +50,8 @@ double computeAnalyticalG(const double x, ...@@ -58,6 +50,8 @@ double computeAnalyticalG(const double x,
return g_analytique; return g_analytique;
} }
// Computes the 'hx' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradHx' should be
// computed. "element" is the element 'j' while will helps to compute 'gradHx'.
double computeAnalyticalGradHx(const double x, double computeAnalyticalGradHx(const double x,
const double y, const double y,
const elementStruct &element) const elementStruct &element)
...@@ -86,6 +80,8 @@ double computeAnalyticalGradHx(const double x, ...@@ -86,6 +80,8 @@ double computeAnalyticalGradHx(const double x,
return gradHx_analytique; return gradHx_analytique;
} }
// Computes the 'hy' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradHy' should be
// computed. "element" is the element 'j' while will helps to compute 'gradHy'.
double computeAnalyticalGradHy(const double x, double computeAnalyticalGradHy(const double x,
const double y, const double y,
const elementStruct &element) const elementStruct &element)
...@@ -114,6 +110,8 @@ double computeAnalyticalGradHy(const double x, ...@@ -114,6 +110,8 @@ double computeAnalyticalGradHy(const double x,
return gradHy_analytique; return gradHy_analytique;
} }
// Computes the 'gx' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradGx' should be
// computed. "element" is the element 'j' while will helps to compute 'gradGx'.
double computeAnalyticalGradGx(const double x, double computeAnalyticalGradGx(const double x,
const double y, const double y,
const elementStruct &element) const elementStruct &element)
...@@ -133,6 +131,8 @@ double computeAnalyticalGradGx(const double x, ...@@ -133,6 +131,8 @@ double computeAnalyticalGradGx(const double x,
return gradGx_analytique; return gradGx_analytique;
} }
// Computes the 'gy' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradGy' should be
// computed. "element" is the element 'j' while will helps to compute 'gradGy'.
double computeAnalyticalGradGy(const double x, double computeAnalyticalGradGy(const double x,
const double y, const double y,
const elementStruct &element) const elementStruct &element)
...@@ -152,7 +152,9 @@ double computeAnalyticalGradGy(const double x, ...@@ -152,7 +152,9 @@ double computeAnalyticalGradGy(const double x,
return gradGy_analytique; return gradGy_analytique;
} }
// Computes the 'g' element numerically. "x" and "y" are the coordinates at which 'g' should be computed. "element" is the
// element 'j' while will helps to compute 'g'. "integrationType" is the type of numerical integration used (in case of numerical
// computation).
double computeNumericalG(const double x, double computeNumericalG(const double x,
const double y, const double y,
const elementStruct &element, const elementStruct &element,
...@@ -182,6 +184,9 @@ double computeNumericalG(const double x, ...@@ -182,6 +184,9 @@ double computeNumericalG(const double x,
return g; return g;
} }
// Computes the 'hx' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradHx' should be computed.
// "element" is the element 'j' while will helps to compute 'gradHx'. "integrationType" is the type of numerical integration used
// (in case of numerical computation).
double computeNumericalGradHx(const double x, double computeNumericalGradHx(const double x,
const double y, const double y,
const elementStruct &element, const elementStruct &element,
...@@ -226,6 +231,9 @@ double computeNumericalGradHx(const double x, ...@@ -226,6 +231,9 @@ double computeNumericalGradHx(const double x,
return gradHx; return gradHx;
} }
// Computes the 'hy' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradHy' should be computed.
// "element" is the element 'j' while will helps to compute 'gradHy'. "integrationType" is the type of numerical integration used
// (in case of numerical computation).
double computeNumericalGradHy(const double x, double computeNumericalGradHy(const double x,
const double y, const double y,
const elementStruct &element, const elementStruct &element,
...@@ -270,6 +278,9 @@ double computeNumericalGradHy(const double x, ...@@ -270,6 +278,9 @@ double computeNumericalGradHy(const double x,
return gradHy; return gradHy;
} }
// Computes the 'gx' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradGx' should be computed.
// "element" is the element 'j' while will helps to compute 'gradGx'. "integrationType" is the type of numerical integration used
// (in case of numerical computation).
double computeNumericalGradGx(const double x, double computeNumericalGradGx(const double x,
const double y, const double y,
const elementStruct &element, const elementStruct &element,
...@@ -299,6 +310,9 @@ double computeNumericalGradGx(const double x, ...@@ -299,6 +310,9 @@ double computeNumericalGradGx(const double x,
return gradGx; return gradGx;
} }
// Computes the 'gy' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradGy' should be computed.
// "element" is the element 'j' while will helps to compute 'gradGy'. "integrationType" is the type of numerical integration used
// (in case of numerical computation).
double computeNumericalGradGy(const double x, double computeNumericalGradGy(const double x,
const double y, const double y,
const elementStruct &element, const elementStruct &element,
......
srcs/BEM/functionsBEM.hpp
+ 199
− 68
View file @ 81ad07f7
This diff is collapsed.
srcs/BEM/mainBEM.cpp
+ 4
− 2
View file @ 81ad07f7
...@@ -11,19 +11,21 @@ int main(int argc, char **argv) ...@@ -11,19 +11,21 @@ int main(int argc, char **argv)
return 0; return 0;
} }
int nthreads = omp_get_max_threads();
omp_set_num_threads(nthreads);
gmsh::initialize(); gmsh::initialize();
gmsh::open(argv[1]); gmsh::open(argv[1]);
gmsh::model::mesh::generate(2); gmsh::model::mesh::generate(2);
// will map the 1d element tag onto the corresponding electroPressure resulting from BEM code
std::map<int, double> electrostaticPressure; std::map<int, double> electrostaticPressure;
int nbViews = 0; int nbViews = 0;
std::map<int,std::pair<double,double>> nodalDisplacements; std::map<int,std::pair<double,double>> nodalDisplacements;
const bool postProcessing = true; const bool postProcessing = true;
const bool meshUntangler = false; const bool meshUntangler = false;
solverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, 1, meshUntangler); // iteration number randomly set to 1 solverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, 1, meshUntangler);
if(postProcessing) if(postProcessing)
gmsh::fltk::run(); gmsh::fltk::run();
......
srcs/BEM/postProcessing.cpp
+ 197
− 133
View file @ 81ad07f7
This diff is collapsed.
srcs/BEM/solverBEM.cpp
+ 37
− 29
View file @ 81ad07f7
#include "functionsBEM.hpp" #include "functionsBEM.hpp"
#include <gmsh.h>
#include <iostream>
#include <map>
#include <sstream>
#include <Eigen/Dense>
#include <cmath>
#include <iomanip>
#include <stdio.h>
#include <chrono>
#include <omp.h>
#include <list>
// Fills the "electrostaticPressure" map that gives the electrostatic pressure associated with each (boundary) element tag. // Fills the "electrostaticPressure" map that gives the electrostatic pressure associated with each (boundary) element tag.
// "nbViews" is the counter of views (useful if other views have to be generated than with this solver). "nodalDisplacements" is // "nbViews" is the counter of views (useful if other views have to be generated than with this solver). "nodalDisplacements" is
...@@ -100,7 +89,9 @@ void solverBEM(std::map<int, double> &electrostaticPressure, ...@@ -100,7 +89,9 @@ void solverBEM(std::map<int, double> &electrostaticPressure,
// Runs the solver for each BEM domain. // Runs the solver for each BEM domain.
for(std::size_t i = 0; i < domainNames.size(); i++) for(std::size_t i = 0; i < domainNames.size(); i++)
singleDomainSolverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, iteration, "BEM_domain_" + domainNames[i], phiViewTag, electricFieldViewTag, computePhi, computeElectricField); singleDomainSolverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, iteration,
"BEM_domain_" + domainNames[i], phiViewTag, electricFieldViewTag, computePhi,
computeElectricField);
if(iteration == 1) if(iteration == 1)
std::cout << "--------------------[BEM SOLVER]--------------------\n\n"; std::cout << "--------------------[BEM SOLVER]--------------------\n\n";
...@@ -135,8 +126,10 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure, ...@@ -135,8 +126,10 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
// Parameters: // Parameters:
bool showTime = false; // true: displays the different time measurements. bool showTime = false; // true: displays the different time measurements.
bool convergence = false; // true: computes the error between the analytical solution and the numerical one on a simple configuration (/!\ use this with the BEM_convergence.geo model). bool convergence = false; // true: computes the error between the analytical solution and the numerical one on a simple
bool visualisation = true; // true: 'ElementNodeData' visualization of the potential, false: 'ElementData' visualization of the potential (Remark: electric field is always displayed with 'ElementData' type). // configuration (/!\ use this with the BEM_convergence.geo model).
bool visualisation = true; // true: 'ElementNodeData' visualization of the potential, false: 'ElementData' visualization of
// the potential (Remark: electric field is always displayed with 'ElementData' type).
bool solutionType = true; // true: analytical computations, false: numerical computations. bool solutionType = true; // true: analytical computations, false: numerical computations.
// Used for displaying the time of different sections. // Used for displaying the time of different sections.
...@@ -152,8 +145,10 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure, ...@@ -152,8 +145,10 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
// Fills the vector of boundary elements. // Fills the vector of boundary elements.
start = omp_get_wtime(); start = omp_get_wtime();
std::map<std::string, std::pair<int, int>> physicalGroups; int domainPhysicalGroupDim;
std::size_t nbDirichletElements = fillElementVector(elementVector, physicalGroups, nodalDisplacements, domainName); int domainPhysicalGroupTag;
std::size_t nbDirichletElements = fillElementVector(elementVector, domainPhysicalGroupDim, domainPhysicalGroupTag,
nodalDisplacements, domainName);
end = omp_get_wtime(); end = omp_get_wtime();
displayTime(start, end, "'fillElementVector'", showTime); displayTime(start, end, "'fillElementVector'", showTime);
...@@ -201,12 +196,12 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure, ...@@ -201,12 +196,12 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
std::vector<double> domainNodeCoords; std::vector<double> domainNodeCoords;
std::vector<int> domainEntityTags; std::vector<int> domainEntityTags;
std::map<int, std::size_t> domainEntityIndices;
std::vector<std::vector<int>> domainElementTypes; std::vector<std::vector<int>> domainElementTypes;
std::vector<std::vector<std::vector<size_t>>> domainElementTags; std::vector<std::vector<std::vector<size_t>>> domainElementTags;
std::map<std::size_t, std::vector<std::size_t>> domainElementNodeTags; std::map<std::size_t, std::vector<std::size_t>> domainElementNodeTags;
int domainPhysicalGroupTag = getModel(domainNodeTags, domainNodeIndices, domainNodeCoords, domainEntityTags, domainEntityIndices, domainElementTypes, domainElementTags, domainElementNodeTags, physicalGroups, domainName); getModel(domainNodeTags, domainNodeIndices, domainNodeCoords, domainEntityTags, domainElementTypes, domainElementTags,
domainElementNodeTags, domainPhysicalGroupDim, domainPhysicalGroupTag);
end = omp_get_wtime(); end = omp_get_wtime();
displayTime(start, end, "'getModel'", showTime); displayTime(start, end, "'getModel'", showTime);
...@@ -224,46 +219,54 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure, ...@@ -224,46 +219,54 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
// Visualisation by ElementNodeData // Visualisation by ElementNodeData
if(visualisation) if(visualisation)
{ {
// Computes a map that stores 'nodeTag' (of the boundary) -> vector of adjacent 'elementIndices' (because only 1 element for higher order elements). // Computes a map that stores 'nodeTag' (of the boundary) -> vector of adjacent 'elementIndices' (because only 1
// element for higher order elements).
std::map<std::size_t, std::vector<std::size_t>> boundaryNodeToElements; std::map<std::size_t, std::vector<std::size_t>> boundaryNodeToElements;
computeNodeToElements(boundaryNodeToElements, elementVector); computeNodeToElements(boundaryNodeToElements, elementVector);
// Computes the value of the potential at each node strictly included in the domain and stores it in the "internalNodalPhi" map // Computes the value of the potential at each node strictly included in the domain and stores it in the
// "internalNodalPhi" map.
start = omp_get_wtime(); start = omp_get_wtime();
std::map<int, double> internalNodalPhi; std::map<int, double> internalNodalPhi;
computeInternalPhi(internalNodalPhi, domainNodeCoords, domainNodeTags, domainNodeIndices, elementVector, integrationType, boundaryNodeToElements, solutionType); computeInternalPhi(internalNodalPhi, domainNodeCoords, domainNodeTags, domainNodeIndices, elementVector,
integrationType, boundaryNodeToElements, solutionType);
end = omp_get_wtime(); end = omp_get_wtime();
displayTime(start, end, "'computeInternalPhi'", showTime); displayTime(start, end, "'computeInternalPhi'", showTime);
// Associates the value of the potential to each node of an element. // Associates the value of the potential to each node of an element.
start = omp_get_wtime(); start = omp_get_wtime();
elementNodeData(elementNodalPhi, domainElementTags, domainElementTypes, domainElementNodeTags, elementVector, internalNodalPhi, domainPhysicalGroupTag, boundaryNodeToElements); elementNodeData(elementNodalPhi, domainElementTags, domainElementTypes, domainElementNodeTags, elementVector,
internalNodalPhi, domainPhysicalGroupTag, boundaryNodeToElements);
end = omp_get_wtime(); end = omp_get_wtime();
displayTime(start, end, "'elementNodeData'", showTime); displayTime(start, end, "'elementNodeData'", showTime);
// Computes the electric field (with "ElementData" type). // Computes the electric field (with "ElementData" type).
start = omp_get_wtime(); start = omp_get_wtime();
computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags, elementVector, solutionType, integrationType, false, computeElectricField); computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags,
elementVector, solutionType, integrationType, false, computeElectricField);
end = omp_get_wtime(); end = omp_get_wtime();
displayTime(start, end, "'computeElementData' (only the electric field)", showTime); displayTime(start, end, "'computeElementData' (only the electric field)", showTime);
// Displays what is needed. // Displays what is needed.
if(computePhi) if(computePhi)
{ {
displayData(phiViewTag, names[0], "results.msh", "ElementNodeData", domainElementTags, elementNodalPhi, viewIndex); displayData(phiViewTag, names[0], "results.msh", "ElementNodeData", domainElementTags, elementNodalPhi,
viewIndex);
viewIndex++; viewIndex++;
} }
if(computeElectricField) if(computeElectricField)
displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, viewIndex); displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags,
elementalElectricField, viewIndex);
} }
// Visualisation by ElementData // Visualisation by ElementData
else else
{ {
// Computes the potential and the electric field. // Computes the potential and the electric field.
start = omp_get_wtime(); start = omp_get_wtime();
computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags, elementVector, solutionType, integrationType, computePhi, computeElectricField); computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags,
elementVector, solutionType, integrationType, computePhi, computeElectricField);
end = omp_get_wtime(); end = omp_get_wtime();
displayTime(start, end, "'computeElementData'", showTime); displayTime(start, end, "'computeElementData'", showTime);
...@@ -274,7 +277,8 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure, ...@@ -274,7 +277,8 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
viewIndex++; viewIndex++;
} }
if(computeElectricField) if(computeElectricField)
displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, viewIndex); displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags,
elementalElectricField, viewIndex);
// Study of the error on a precise model. // Study of the error on a precise model.
if(convergence) if(convergence)
...@@ -282,9 +286,13 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure, ...@@ -282,9 +286,13 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
std::cout << "\tCONVERGENCE:\n"; std::cout << "\tCONVERGENCE:\n";
std::string convergenceModelName = "BEM_convergence"; std::string convergenceModelName = "BEM_convergence";
if(names[0].compare(convergenceModelName) == 0) if(names[0].compare(convergenceModelName) == 0)
convergenceFun(domainEntityTags, domainElementTypes, domainElementTags, elementalPhi, elementalElectricField); convergenceFun(domainEntityTags, domainElementTypes, domainElementTags, elementalPhi,
elementalElectricField);
else else
std::cout << "\tWarning: Cannot check convergence on model '" << names[0] << ".geo'. Use the '" << convergenceModelName + ".geo" << "' model if you want to check for convergence.\n\n"; {
std::cout << "\tWarning: Cannot check convergence on model '" << names[0] << ".geo'. Use the '";
std::cout << convergenceModelName + ".geo" << "' model if you want to check for convergence.\n\n";
}
} }
} }
} }
......
srcs/BEM/verification.cpp
+ 11
− 35
View file @ 81ad07f7
#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 "functionsBEM.hpp"
#include <omp.h>
#include <algorithm>
#include <list>
// Displays the time taken between the values of "start" and "end" (in seconds) if "showTime" is set to 'true'. "functionName" is
// the name of the function that has been timed.
void displayTime(const double start, void displayTime(const double start,
const double end, const double end,
std::string functionName, const bool showTime) std::string functionName, const bool showTime)
{ {
//displays the time of the desired code section
if(showTime) if(showTime)
{ {
double duration = end - start; double duration = end - start;
std::cout << functionName << "\n"; std::cout << functionName << "\n";
std::cout << "Elapsed time: " << 1000000*duration << " [µs]\tor " << duration << " [s]\n\n"; std::cout << "Elapsed time: " << 1000000*duration << " [µs]\tor " << duration << " [s]\n\n";
// std::cout << std::endl;
// std::cout << function_name << std::endl;
// std::cout << "Elapsed time in nanoseconds: "
// << std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count()
// << " ns" << std::endl;
// std::cout << "Elapsed time in microseconds: "
// << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
// << " µs" << std::endl;
// std::cout << "Elapsed time in milliseconds: "
// << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count()
// << " ms" << std::endl;
// std::cout << "Elapsed time in seconds: "
// << std::chrono::duration_cast<std::chrono::seconds>(end - start).count()
// << " sec";
// std::cout << std::endl;
// std::cout << std::endl;
} }
} }
// Displays the different errors induced by the boundary element method compared with a simple analytical solution on a certain
// model geometry ("BEM_convergence.geo"). "domainEntityTags" is a vector containing the tag of each entity of the domain.
// "domainElementTypes" stores the types of the elements of the the domain. The first dimension of this vector represents the
// different entities and the second, the different element types. "domainElementTags" stores the tags of the elements of the
// domain. The first dimension of this vector represents the different entities, the second, the different element types and the
// third, the different element tags. "elementalPhi" and "elementalElectricField" are the used data for comparing numerical and
// analytical solutions. The convergence function can only be used with the "ElementData" data type.
void convergenceFun(const std::vector<int> &domainEntityTags, void convergenceFun(const std::vector<int> &domainEntityTags,
const std::vector<std::vector<int>> &domainElementTypes, const std::vector<std::vector<int>> &domainElementTypes,
const std::vector<std::vector<std::vector<size_t>>> &domainElementTags, const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
......
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