Skip to content
Snippets Groups Projects
Commit a8b3721e authored by Denis Louis's avatar Denis Louis
Browse files

Kevin update

parent 65678854
No related branches found
No related tags found
1 merge request!8SECOND PROGRESS DEADLINE
This diff is collapsed.
#include <gmsh.h>
#include <iostream>
#include <vector> // Why should I include this ?
#include <map>
#include <sstream>
#include <Eigen/Dense>
#include <cmath>
#include <iomanip>
#include <stdio.h>
#include <chrono>
#include <omp.h>
// Each element has all its information gathered as in this structure.
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]
std::size_t nodeTags [2]; // tags of the [first, second] nodes of the element
size_t nodeTags [2]; // tags of the [first, second] nodes of the element
size_t midNodeTags;
double length;
double x1, y1, x2, y2;
double midX, midY;
double deltaX, deltaY;
};
int solverBEM(int argc, char **argv, std::map<int, double> &electrostaticPressure, int & nbViews);
void dispNodeCoordFun(const std::vector<size_t> &allNodeTags, const std::vector<double> &allCoord);
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(std::map<size_t, int> allNodeIndex, const std::vector<double> &allCoord, std::vector<elementStruct> &elementVector, int &dirichletCount, const bool dispHierarchy, std::vector<size_t> &NodeTags2d_domain, std::vector<double> &NodeCoord2d_domain, std::vector<int> &elementTypes2d_domain_BEM, std::vector<std::vector<size_t>> &elementTags2d_domain_BEM, std::vector<std::vector<size_t>> &elnodeTags2d_domain_BEM, int &entities_tags_for_physical_group_domain_BEM, std::vector<int> &entities_tags_2d_domain_BEM);
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);
void computeData(const std::vector<elementStruct> &elementVector, std::vector<size_t> &tags, std::vector<std::vector<double>> &data, const std::string integrationType, std::vector<int> &entities_tags_2d_domain_BEM);
void displayData(const std::string viewName, const std::string modelName, const std::string fileName, const std::vector<size_t> &tags, const std::vector<std::vector<double>> &data);
void compute_internal_phi(std::vector<double> &allCoord,std::vector<size_t> &NodeTags2d_domain, std::vector<elementStruct> &elementVector, std::map<int,double> &nodeMapDomain, const std::string integrationType);
void element_node_data(std::vector<std::vector<std::vector<double>>> &data, std::vector<std::vector<size_t>> &elementTags2d_domain_BEM, std::vector<int> &elementTypes2d_domain_BEM, std::vector<std::vector<size_t>> &nodeTags2d_domain_BEM, std::vector<elementStruct> &elementVector, std::map<int, double> &nodeMapDomain, bool &show_element_node_data, int &entities_tags_for_physical_group_domain_BEM);
void disp_element_node_data(std::vector<int> &nodes_commun_tag, std::vector<int> &elements_commun_tag);
void display_time(std::chrono::time_point<std::chrono::system_clock> &start, std::chrono::time_point<std::chrono::system_clock> &end, std::string &function_name, bool &show_time);
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);
void computeData(const std::vector<elementStruct> &elementVector, const std::vector<int> &elementTypes, const std::vector<std::vector<size_t>> &elementTags, std::vector<size_t> &tags, std::vector<std::vector<double>> &data);
void displayData(const std::string viewName, const std::string modelName, const std::string fileName, const std::vector<size_t> &tags, const std::vector<std::vector<double>> &data);
int solverBEM(int argc, char **argv, std::map<int, double> &electrostaticPressure, int &nbViews);
\ No newline at end of file
......@@ -26,13 +26,13 @@ int main(int argc, char **argv)
std::map<int, double>::iterator it;
for (it = electrostaticPressure.begin(); it != electrostaticPressure.end(); it++)
/*for (it = electrostaticPressure.begin(); it != electrostaticPressure.end(); it++)
{
std::cout << it->first // string (key)
<< ':'
<< it->second // string's value
<< std::endl;
}
}*/
gmsh::fltk::run();
......
#include "functionsBEM.hpp"
#include <gmsh.h>
#include <iostream>
#include <map>
#include <sstream>
#include <math.h>
#include <Eigen/Dense>
#include <cmath>
#include <iomanip>
#include <stdio.h>
#include <chrono>
#include <omp.h>
#include <list>
int solverBEM(int argc, char **argv, std::map<int, double> &electrostaticPressure, int & nbViews)
{
......@@ -16,29 +20,71 @@ int solverBEM(int argc, char **argv, std::map<int, double> &electrostaticPressur
}
// Parameters for displaying variables (for debugging).
bool dispNodeCoord = false; // Coordinates of the nodes.
bool dispHierarchy = false; // Hierarchy (physical groups -> entities -> elements -> nodes).
bool dispMatrices = false; // Matrices A, B, c, b, x. /!\ May lead to visualisation problems for large matrices.
bool dispValues = false; // Elements and their values of u and u_n.
bool openGUI = false; // Graphical interface of GMSH.
bool dispNodeCoord = 0; // Coordinates of the nodes.
bool dispHierarchy = 0; // Hierarchy (physical groups -> entities -> elements -> nodes).
bool dispMatrices = 0; // Matrices A, B, c, b, x. /!\ May lead to visualisation problems for large matrices.
bool dispValues = 0; // Elements and their values of u and u_n.
bool openGUI = 0; // Graphical interface of GMSH.
bool show_element_node_data = 0;
bool show_time = 0; //display the elapsed time for the internal potential calculus in the terminal
int visualisation = 1; //0: no visualisation, 1: element node data, 2: element data
//to display the time of the different functions
std::chrono::time_point<std::chrono::system_clock> start, end;
std::chrono::time_point<std::chrono::system_clock> start_total, end_total;
std::string function_name;
start_total = std::chrono::system_clock::now();
// Initialises a vector of elements (that will contain all the elements of interrest)
std::vector<elementStruct> elementVector;
std::string integrationType = "Gauss4"; // Could vary depending on the distance between the elements.
// 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::map<size_t, int> allNodeIndex; // nodeTag -> position in allNodeTags and allCoord
std::vector<size_t> allNodeTags;
std::vector<double> allCoord;
std::vector<double> allParametricCoord;
gmsh::model::mesh::getNodes(allNodeTags, allCoord, allParametricCoord);
for(size_t i = 0; i < allNodeTags.size(); i++)
{
allNodeIndex[allNodeTags[i]] = i;
}
// Print the coordinates of all nodes if dispNodeCoord != 0
if(dispNodeCoord)
{
dispNodeCoordFun(allNodeTags, allCoord);
}
// Initially, these are 0 and are incremented each time a dirichlet/neumann-element is encountered.
int dirichletCount = 0;
getBC(allNodeIndex, allCoord, elementVector, dirichletCount, dispHierarchy);
//pour compute_internal_phi
std::vector<size_t> NodeTags2d_domain;
std::vector<double> NodeCoord2d_domain;
// vector<double> NodeParametricCoord2d_domain;
// model::mesh::getNodes(NodeTags2d_domain, NodeCoord2d_domain, NodeParametricCoord2d_domain,2);
//pour elementnodedata
std::vector<int> elementTypes2d_domain_BEM;
std::vector<std::vector<size_t> > elementTags2d_domain_BEM;
std::vector<std::vector<size_t> > nodeTags2d_domain_BEM;
int entities_tags_for_physical_group_domain_BEM;
std::vector<int> entities_tags_2d_domain_BEM;
// model::mesh::getElements(elementTypes2d,elementTags2d,nodeTags2d,2);
start = std::chrono::system_clock::now();
//calculates the potential at the barycentre of each element of the domain
getBC(allNodeIndex,allCoord, elementVector, dirichletCount, dispHierarchy, NodeTags2d_domain, NodeCoord2d_domain,elementTypes2d_domain_BEM,elementTags2d_domain_BEM,nodeTags2d_domain_BEM,entities_tags_for_physical_group_domain_BEM, entities_tags_2d_domain_BEM);
end = std::chrono::system_clock::now();
function_name = "getBC";
display_time(start, end, function_name, show_time);
std::vector<std::string> names;
gmsh::model::list(names);
//std::cout << "the file contains " << names.size() << " model names\n";
// 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());
......@@ -46,12 +92,17 @@ int solverBEM(int argc, char **argv, std::map<int, double> &electrostaticPressur
Eigen::MatrixXd c = Eigen::MatrixXd::Zero(elementVector.size(), 1);
// Double loop over each element in order to compute the element h and g that will be given to the matrices A and B.
start = std::chrono::system_clock::now();
double h;
double g;
#pragma omp parallel for schedule(guided) shared(elementVector,A,B,c,integrationType,dirichletCount) private(h,g)
for(int i = 0; i < (int) elementVector.size(); i++)
{
for(int j = 0; j < (int) elementVector.size(); j++)
{
// Will be computed and stored in A and B.
double h;
double g;
h = 0;
g = 0;
if(i == j) // Formula for diagonal elements.
{
......@@ -61,9 +112,9 @@ int solverBEM(int argc, char **argv, std::map<int, double> &electrostaticPressur
else // General case
{
h = computeH(elementVector[i].midX, elementVector[i].midY, elementVector[j]);
g = computeG(elementVector[i].midX, elementVector[i].midY, elementVector[j], "Gauss4");
g = computeG(elementVector[i].midX, elementVector[i].midY, elementVector[j], integrationType);
}
// Storage of h and g at the right place (in matrices A and B).
if(j < dirichletCount)
{
......@@ -78,83 +129,124 @@ int solverBEM(int argc, char **argv, std::map<int, double> &electrostaticPressur
c(j) = elementVector[j].bcValue[1];
}
}
}
end = std::chrono::system_clock::now();
function_name = "compute H and G";
display_time(start, end, function_name, show_time);
//If we want to display the time it takes for eigen to solve the linear system
start = std::chrono::system_clock::now();
Eigen::MatrixXd b = B*c;
Eigen::MatrixXd x = A.fullPivLu().solve(b);
end = std::chrono::system_clock::now();
function_name = "solve linear system";
display_time(start, end, function_name, show_time);
// Matrices verification
if(dispMatrices)
{
dispMatricesFun(A, B, c, b, x);
}
// Assign each value of x to the right element.
for(int i = 0; i < (int) elementVector.size(); i++)
{
if(i < dirichletCount)
if( i < dirichletCount)
elementVector[i].bcValue[1] = x(i);
else
elementVector[i].bcValue[0] = x(i);
}
if(dispValues)
{
dispValuesFun(elementVector);
}
/*
* Visualisation
*/
std::vector<std::string> names;
gmsh::model::list(names);
//std::vector<std::vector<double>> phi(nb2dElements);
//std::vector<std::size_t> tags(nb2dElements);
std::vector<std::vector<double>> phi;
std::vector<std::size_t> tags;
// get physical groups
std::map<std::string, std::pair<int, int>> groups;
gmsh::vectorpair dimTags;
gmsh::model::getPhysicalGroups(dimTags);
for (std::size_t i = 0; i < dimTags.size(); ++i)
//visualisation by elementnodedata
if(visualisation == 1)
{
int dim = dimTags[i].first;
int tag = dimTags[i].second;
std::string name;
gmsh::model::getPhysicalName(dim, tag, name);
groups[name] = {dim, tag};
}
// get all the 2D elements in the BEM domain
std::string groupname = "BEM_domain";
int dim = groups[groupname].first;
int tag = groups[groupname].second;
std::vector<int> entities_tags;
gmsh::model::getEntitiesForPhysicalGroup(dim, tag, entities_tags);
std::vector<int> elementTypes;
std::vector<std::vector<std::size_t>> elementTags;
std::vector<std::vector<std::size_t>> elnodeTags;
int nb2dElements = 0;
for(std::size_t i=0; i< entities_tags.size(); ++i)
{
gmsh::model::mesh::getElements(elementTypes, elementTags, elnodeTags, dim, entities_tags[i]);
for(std::size_t j=0; j< elementTags.size(); ++j)
//The first index indicates the elementType, the second the elementTag and the third the nodesTags associated with this element. For example, it says whether it is a square or a triangle, the element tag and finally the tag of the 3 or 4 nodes associated with it.
std::vector<std::vector<std::vector<double>>> data;
//This map associates the value of the potential to all nodes strictly included in the domain and not on the boundary.
std::map<int, double> nodeMapDomain;
//To displays the time of the this code section
//the purpose of this function is to fill the nodemapdomain map and thus to calculate the value of the potential at each node in the domain
start = std::chrono::system_clock::now();
//calculates the potential at the barycentre of each element of the domain
compute_internal_phi(allCoord,NodeTags2d_domain,elementVector,nodeMapDomain,integrationType);
end = std::chrono::system_clock::now();
function_name = "compute_internal_phi";
display_time(start, end, function_name, show_time);
//this function allows to correctly associate the value of the potential to each node of an element in order to display the ElementNodeData view
start = std::chrono::system_clock::now();
element_node_data(data,elementTags2d_domain_BEM,elementTypes2d_domain_BEM,nodeTags2d_domain_BEM,elementVector, nodeMapDomain,show_element_node_data,entities_tags_for_physical_group_domain_BEM);
end = std::chrono::system_clock::now();
std::string function_name = "ElementNodeData";
display_time(start, end, function_name, show_time);
int viewNodeTag = gmsh::view::add("phi");
for(size_t i = 0; i < data.size(); i++)
{
nb2dElements += elementTags[j].size();
gmsh::view::addModelData(viewNodeTag,0,names[0],"ElementNodeData",elementTags2d_domain_BEM[i],data[i],0.0,1);
}
computeData(elementVector, elementTypes, elementTags, tags, phi);
}
gmsh::view::write(viewNodeTag, "results.msh");
displayData("phi", names[0], "results.msh", tags, phi);
gmsh::option::setNumber("View[0].IntervalsType", 3);
//option::setNumber("View[0].AdaptVisualizationGrid", 1);
gmsh::option::setNumber("View[0].MaxRecursionLevel", 3);
gmsh::option::setNumber("View[0].TargetError", -0.0001);
}
//visualisation by element barycenter
if(visualisation == 2)
{
std::vector<std::string> names;
gmsh::model::list(names);
//To get the total number of 2d elements that make up the domain, we simply add the number of squares to the number of triangles,...
int nb2dElements = 0;
for(size_t i = 0; i < elementTypes2d_domain_BEM.size(); i++)
nb2dElements += elementTags2d_domain_BEM[i].size();
//Tag and associated potential of all 2d elements. In this view the potential is calculated at the barycentre and is constant per element.
std::vector<std::vector<double>> phi(nb2dElements,std::vector<double>(1,0));
std::vector<size_t> tags(nb2dElements);
//To displays the time of the this code section
start = std::chrono::system_clock::now();
//calculates the potential at the barycentre of each element of the domain
computeData(elementVector, tags, phi, integrationType, entities_tags_2d_domain_BEM);
end = std::chrono::system_clock::now();
function_name = "computeData";
display_time(start, end, function_name, show_time);
//For ease of use, a function handles the code for the display.
displayData("phi", names[0], "results.msh", tags, phi);
}
//Calculates the electrostatic pressure per element in order to provide these results to the FEM code that will handle the rest of the algorithm
double epsilon = 8.8541878128e-12;
for(size_t i = 0; i < elementVector.size(); i++)
{
electrostaticPressure[elementVector[i].tag] = epsilon * pow(elementVector[i].bcValue[1], 2.0);
}
end_total = std::chrono::system_clock::now();
function_name = "total time of the code BEM";
display_time(start_total, end_total, function_name, show_time);
if(openGUI)
{
gmsh::fltk::run();
}
nbViews = 1; // number of views generated (useful for FEM display of results)
return 0;
}
\ No newline at end of file
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment