Skip to content

Commit

Permalink
Nothing works...
Browse files Browse the repository at this point in the history 
This code architecture is not suitable:

- Split into classes reduce performances and does not
  improve readability.
- Eigen library is a nightmare for large scale project.
- Does not take advantages of GMSH api, cache, ...
  About everything is stored 2 times in memory. (no kidding)
- Debug nightmares and Spaghetti coding (class dependencies).
  • Loading branch information
@povanberg
povanberg committed Mar 4, 2019
1 parent 1bc5758 commit d68cf77
Show file tree
Hide file tree
Showing 7 changed files with 310 additions and 43 deletions.
11 changes: 8 additions & 3 deletions include/Element.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ class Element {
std::vector<Eigen::VectorXd> basisFcts; // [f](g) : f= basis fct; g=int point
std::vector<Eigen::MatrixXd> gradBasisFcts; // [f](g,i) : f= fct; g=int point; i=df/dx_i
//----------------------------------------------------------------------------------------------------------
Eigen::VectorXd u;
public:
Element(int dim, int tag, std::vector<int> &nodeTags);
std::vector<Face> faces;
Expand All @@ -38,6 +39,7 @@ class Element {
const int &getType();
const int &getOrder();
const int &getNumNodes();
std::vector<int> &getNodeTags();
const std::string &getName();
const Eigen::Matrix3d &getJacobian(const int g);
const Eigen::Matrix3d &getInvJacobian(const int g);
Expand All @@ -48,12 +50,15 @@ class Element {
const Eigen::VectorXd &getBasisFcts(const int f);
const Eigen::MatrixXd &getGradBasisFcts(const int f);
void getMassMatrix(Eigen::MatrixXd &massMAtrix);
void getStiffMatrix(Eigen::MatrixXd &stiffMatrix, const Eigen::Vector3d a);
void getStiffMatrix(Eigen::MatrixXd &stiffMatrix, const Eigen::Vector3d &a);
void getData(Eigen::VectorXd &data);
void setData(Eigen::VectorXd &data);
//----------------------------------------------------------------------------------------------------------
// Check if element contains the face
// with corresponding tag

bool hasNode(const int tag);

void getFlux(Eigen::MatrixXd &Flux, const Eigen::Vector3d &a);

// Add face to element
Element &addFace(Face face);

Expand Down
31 changes: 26 additions & 5 deletions include/Face.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,27 +7,48 @@

class Face {
private:
//----------------------------------------------------------------------------------------------------------
int tag; // Tag of the face
int dim; // Dimension of the face
int numNodes; // Number of nodes on surface
int type; // Type of the face
int numIntPoints; // Number of integration points
int numBasisFcts; // Number of basis fcts
std::string name; // Face Type name
std::vector<int> nodeTags; // List of nodes of the face
std::vector<int> elementTags; // List of adjacent element tags
Eigen::Vector3d normal; // Normal (i) : n_i
//----------------------------------------------------------------------------------------------------------
void setNormal(); // Compute normal of face
Eigen::MatrixXd basisFcts; // (f,g) : f=basis fct, g=gauss
Eigen::MatrixXd xPoints; // (g,i) : g=int point; x_i
Eigen::MatrixXd uPoints; // (g,i) : g=int point; u_i
Eigen::VectorXd weights; // (g) : g=int point
std::vector<Eigen::Matrix3d> jacobian; // [g](i,j) : g=int point; dx_i/du_j
Eigen::VectorXd detJacobian; // (g) : g=int point

public:
Face(int tag, std::string name, int dim, int numNodes, int type, std::vector<int> &nodeTags);
//----------------------------------------------------------------------------------------------------------

bool boundary;

void setNormal();
void setBasisFcts();
const int &getTag();
const std::vector<int> &getNodeTags();
const Eigen::Vector3d &getNormal();
//----------------------------------------------------------------------------------------------------------

bool hasNode(const int tag);

double getFluxInt(const Eigen::Vector3d &a);
int getSecondElement(const int tag1);

// Add adjacent element to the face
Face &addElement(int tag);

// Set Jacobian
void setJacobian(std::vector<double> &jacobian,
std::vector<double> &detJacobian,
std::vector<double> &xPoints);


};


Expand Down
16 changes: 12 additions & 4 deletions include/Mesh.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,21 +9,29 @@
class Mesh {

private:
//----------------------------------------------------------------------------------------------------------

int dim; // Mesh dimension
std::string name; // File name
std::map<std::string, std::pair<int,int>> physDimTags; // Associate physical name to tag/dim
std::vector<int> domainEntityTags; // List of all geometrical entities in domain
std::vector<int> diricheletEntityTags; // List of all geometrical entities in dirichelet
//----------------------------------------------------------------------------------------------------------

public:

Mesh(std:: string name);

std::vector<Element> elements; // List of elements in Mesh
Element &getElement(const int tag);

void getMassMatrix(Eigen::SparseMatrix<double> &M);
void getStiffMatrix(Eigen::SparseMatrix<double> &K, const Eigen::Vector3d &a);
void getFlux(Eigen::VectorXd &F, const Eigen::Vector3d &a);

Element &getDataElement(const int tag);
void getNodeTags(std::vector<int> &nodeTags);
void getData(Eigen::VectorXd &data);
void setData(Eigen::VectorXd &data);

void getMassMatrix(Eigen::SparseMatrix<double> M);
void getStiffMatrix(Eigen::SparseMatrix<double> K, const Eigen::Vector3d a);
};


Expand Down
33 changes: 32 additions & 1 deletion src/Element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ Element::Element(int dim, int tag, std::vector<int> &nodeTags) {
this->order,
this->numNodes,
this->paramCoord);
this->u.resize(this->numNodes);
}

// Add face to the element
Expand Down Expand Up @@ -119,6 +120,11 @@ const std::string &Element::getName(){
return this->name;
}

// Get list of node tags
std::vector<int> &Element::getNodeTags(){
return this->nodeTags;
}

// Get jacobian at gauss node g
const Eigen::Matrix3d &Element::getJacobian(const int g){
return this->jacobian[g];
Expand Down Expand Up @@ -177,7 +183,7 @@ void Element::getMassMatrix(Eigen::MatrixXd &massMatrix){
}

// Get element mass matrix
void Element::getStiffMatrix(Eigen::MatrixXd &stiffMatrix, const Eigen::Vector3d a){
void Element::getStiffMatrix(Eigen::MatrixXd &stiffMatrix, const Eigen::Vector3d &a){
stiffMatrix.resize(this->numBasisFcts, this->numBasisFcts);
for(unsigned int i=0; i<this->numBasisFcts; ++i){
for(unsigned int j=0; j<this->numBasisFcts; ++j){
Expand All @@ -186,4 +192,29 @@ void Element::getStiffMatrix(Eigen::MatrixXd &stiffMatrix, const Eigen::Vector3d
.dot(detJacobian);
}
}
}

void Element::getFlux(Eigen::MatrixXd &Flux, const Eigen::Vector3d &a){
// Independent of numerical flux
Flux = Eigen::MatrixXd::Zero(this->numBasisFcts, this->numBasisFcts);
for(unsigned int i=0; i<this->numBasisFcts; ++i) {
for (unsigned int j = 0; j < this->numBasisFcts; ++j) {
for(Face &face: this->faces){
if(face.hasNode(this->nodeTags[i]) && face.hasNode(this->nodeTags[j])){
Flux(i, j) += face.getFluxInt(a);
}
else{
Flux(i, j) += 0.0;
}
}
}
}
}

void Element::getData(Eigen::VectorXd &data){
data = this->u;
}

void Element::setData(Eigen::VectorXd &data){
this->u = data;
}
63 changes: 63 additions & 0 deletions src/Face.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ Face::Face(int tag, std::string name, int dim, int numNodes, int type, std::vect
this->type = type;
this->nodeTags = nodeTags;
setNormal();
setBasisFcts();
}

void Face::setNormal(){
Expand All @@ -35,12 +36,74 @@ void Face::setNormal(){
}
}

void Face::setBasisFcts() {
// Gmsh api call
int numComp;
std::vector<double> basisFcts;
std::vector<double> uPoints;
gmsh::model::mesh::getBasisFunctions(this->type, "Gauss3", "Lagrange",
uPoints, numComp, basisFcts);
// Eigen conversion
this->numBasisFcts = this->numNodes;
this->numIntPoints = (int) basisFcts.size() / this->numBasisFcts;
this->basisFcts = Eigen::Map<Eigen::MatrixXd>(basisFcts.data(), this->numBasisFcts, this->numIntPoints);

this->weights = Eigen::VectorXd(this->numIntPoints);
this->uPoints = Eigen::MatrixXd(this->numIntPoints, 3);
for(int g = 0; g < this->numIntPoints; ++g){
this->uPoints(g, 0) = uPoints[g*4+0];
this->uPoints(g, 1) = uPoints[g*4+1];
this->uPoints(g, 2) = uPoints[g*4+2];
this->weights(g) = uPoints[g*4+3];
}
}

void Face::setJacobian(std::vector<double> &jacobian,
std::vector<double> &detJacobian,
std::vector<double> &xPoints){

this->numIntPoints = numIntPoints;
this->numBasisFcts = this->numNodes;
std::vector<double>::const_iterator gIt = jacobian.begin();
for(unsigned int g=0; g<this->numIntPoints; ++g, gIt+=9){
std::vector<double> gJacobian(gIt, gIt + 9);
this->jacobian.push_back(Eigen::Map<Eigen::Matrix3d>(gJacobian.data()).transpose());
}
this->detJacobian = Eigen::Map<Eigen::VectorXd>(detJacobian.data(), this->numIntPoints);
this->xPoints = Eigen::Map<Eigen::MatrixXd>(xPoints.data(), 3, this->numIntPoints).transpose();
}

// Add adjacent to the face
Face &Face::addElement(int tag) {
this->elementTags.push_back(tag);
if(this->elementTags.size() > 1)
this->boundary = false;
else
this->boundary = true;
return *this;
}

bool Face::hasNode(const int tag){
auto it = std::find(this->nodeTags.begin(), this->nodeTags.end(), tag);
if(it != this->nodeTags.end())
return true;
else
return false;
}

int Face::getSecondElement(const int tag1){
if(tag1 == this->elementTags[0])
return this->elementTags[1];
else
return this->elementTags[0];
}

double Face::getFluxInt(const Eigen::Vector3d &a){
return a.dot(this->normal)*weights.cwiseProduct(basisFcts.row(0).transpose())
.cwiseProduct(basisFcts.row(1).transpose())
.dot(detJacobian);
}

// Get Nodes tags
const std::vector<int> &Face::getNodeTags(){
return this->nodeTags;
Expand Down
Loading

0 comments on commit d68cf77

Please sign in to comment.