FNELEM-GPU

Simple finite elements matrix structural analysis using GPU. Actually supports 2D planar membrane/shell.

Requirements

C++ STANDARD 14 and CUDA. No other libraries are needed.

Import

To fully use this library, you must include the following files:

#include "fnelem/math/fematrix.cpp"
#include "fnelem/math/fematrix_utils.cpp"
#include "fnelem/math/matrix_inversion_cpu.cpp"
#include "fnelem/math/matrix_inversion_cuda.cu"

#include "fnelem/analysis/static_analysis.cpp"
#include "fnelem/model/base/model.cpp"
#include "fnelem/model/base/model_component.cpp"
#include "fnelem/model/elements/element.cpp"
#include "fnelem/model/elements/membrane.cpp"
#include "fnelem/model/loads/load.cpp"
#include "fnelem/model/loads/load_membrane_distributed.cpp"
#include "fnelem/model/loads/load_node.cpp"
#include "fnelem/model/loads/load_pattern.cpp"
#include "fnelem/model/loads/load_pattern_constant.cpp"
#include "fnelem/model/nodes/node.cpp"
#include "fnelem/model/restraints/restraint.cpp"
#include "fnelem/model/restraints/restraint_node.cpp"

Usage

Model

Model class is the base object for defining a problem. It needs the number of degrees of freedom, the dimension (actually only 2D is supported), the nodes, elements, restraints, loads and load patterns.

Model *model = new Model(dimension, number_degrees_of_freedom);

// Define nodes
model->add_nodes(nodes);
model->add_elements(elements);
model->add_restraints(restraints);
model->add_load_patterns(load_pattern);

Nodes, elements, restraints and load_pattern are std::vector objects that contains the library Objects.

StaticAnalysis

StaticAnalysis class performs the basic static analysis method, using the model definition, it creates the structure stiffness matrix, then using K*u = Fit solves the displacement vector of the nodes.

StaticAnalysis *analysis = new StaticAnalysis(model);
analysis->analyze(use_gpu);

If use_gpu is true then GPU performs the main matrix inversion. After the analysis is complete, the model can save the results into a file. That file contains the main structure elements: nodes, shells, reactions and internal forces of the elements.

model->save_results("file.txt");

Example of the output:

FNELEM-GPU -  Finite element structural analysis using CUDA and GPU.
              v0.2.0 (21/12/2018) @ ppizarror

--------------------------------------------------------------------
Input model properties:
--------------------------------------------------------------------

Nodes:
	Node count:	6
	Node N1:	0.000000	0.000000
	...

Elements:
	Element count:		2
	Membrane MEM1:
		Width (2b):		2
		Height (2h):	2
		Thickness:		1
		Elastic mod:	2000
		Poisson mod:	0.2
		Element nodes:	N1, N3, N4, N2
	...

--------------------------------------------------------------------
Analysis results:
--------------------------------------------------------------------

Node displacements:
	Node N1:	0.000000	0.000000
	...

Node reactions:
	Node N1:	369.158879	500.000000
	...

Element stresses:
	Membrane MEM1:
		Node N1 (-b, -h):	369.159,	500
		Node N3 (+b, -h):	130.841,	0
		Node N4 (+b, +h):	-500,	-500
		Node N2 (-b, +h):	0,	0
		Stress MEM1 [GLX GLY X Y SIGMAX SIGMAY SIGMAXY DISPLX DISPLY]
			0 0 -1 -1 -44.3925 -221.963 -193.9250 0
			...

Node

The node class must be defined using the global coordinates. Nodes are used by the elements.

std::vector<Node *> *nodes = new std::vector<Node *>();
nodes->push_back(new Node(node_tag, coordinate_x, coorinate_y));
nodes->push_back(new Node(node_tag, coordinate_x, coordinate_y, coordinate_z));

Element

Membrane shell must be defined using four nodes, the elastic and poisson modulus and the thickness.

std::vector<Element *> *elements = new std::vector<Element *>();
elements->push_back(new Membrane(membrane_tag, node1, node2, node3, node4, E, poisson, thickness));

Restraints

The restraint class must be used when a node cannot move in any direction.

std::vector<Restraint *> *restraints = new std::vector<Restraint *>();
RestraintNode *r = new RestraintNode(restraint_tag, node);
r1->add_dofid(index);
r1->add_all(); // This method adds all degrees of freedom
restraints->push_back(r);

Load

LoadNode

NodeLoad is a load applied into a node, it needs the node reference and a load vector. In FNELEM-GPU matrices and vectors are defined using FEMatrix class. That class provides the following main methods:

• Add
• Substract
• Multiply
• Inverse (Using CPU and GPU)
• Determinant
• max(), min()
• clone()

To define a vector:

FEMatrix *vector = FEMatrix_vector(size);
vector->set(index, value);
vector->get(index); 

To define a matrix:

FEMatrix *matrix = FEMatrix_vector(n, m);
matrix->set(i, j, value);
matrix->get(i, j); 

Then, node load is defined by:

// Create load vector
std::vector<Load *> *loads = new std::vector<Load *>();

// Create LoadNode
FEMatrix *load = FEMatrix_vector(size);
load->set(index, value);
...
loads->push_back(new LoadNode(load_tag, node, load));

LoadMembraneDistributed

This class can apply a distributed load. It needs the membrane, the node position (1, 2, 3, 4)

std::vector<Load *> *loads = new std::vector<Load *>();
loads->push_back(new LoadMembraneDistributed(load_tag, membrane, position_initial,
    position_final, load_at_initial, distance_from_initial, load_at_final, distance_from_initial));

Load pattern

The load pattern must be used to apply the loads into the model. Actually only the constant load pattern is defined.

std::vector<LoadPattern *> *loadpattern = new std::vector<LoadPattern *>();
loadpattern->push_back(new LoadPatternConstant(load_tag, loads));

Example

The following code defines a simple bridge, composed by N membranes, under a distributed load of 100kN/m. The membranes are the same, 100cm height, 100cm width, 15cm thickness, elastic modulus of 300000kN/cm^2, Poisson modulus of 0.2. This example can also be found in test/analysis/test_static_analysis.h.

int N = 3; // Number of membranes
int b = 100; // Pier width
int h = 100; // Pier height

//    N = 2
//       4 ------ 5 ------ 6
//       |        |        |
//     h |   (1)  |   (2)  |  ......
//       |        |        |
//       1 ------ 2 ------ 3 .....
//       ^    b            ^
//
//    N = 3
//       4 ------ 5 ------ 6 ------- 8
//       |        |        |         |
//     h |   (1)  |   (2)  |   (3)   |
//       |        |        |         |
//       1 ------ 2 ------ 3 ------- 4
//       ^    b                      ^
// =====================================

double t = 15; // Thickness (cm)
double E = 300000; // Elastic modulus
double nu = 0.15; // Poisson modulus

// Number degrees of freedom
int gdl = N * 4;

// Create model
Model *model = new Model(2, gdl);

// Create nodes
std::vector<Node *> *nodes = new std::vector<Node *>();
int j;
for (int i = 1; i <= N + 1; i++) {
    nodes->push_back(new Node("N" + std::to_string(i), b * (i - 1), 0));
}
for (int i = 1; i <= N + 1; i++) {
    j = N + 1 + i;
    nodes->push_back(new Node("N" + std::to_string(j), b * (i - 1), h));
}

// Add nodes to model
model->add_nodes(nodes);

// Create elements
// n4 ------------ n3
//  |              |
//  |      (i)     |
//  |              |
// n1 ------------ n2
unsigned long n1, n2, n3, n4;
std::vector<Element *> *elements = new std::vector<Element *>();
for (unsigned long i = 0; i < N; i++) {
    n1 = i;
    n2 = i + 1;
    n3 = N + i + 2;
    n4 = N + i + 1;
    elements->push_back(new Membrane("MEM" + std::to_string(i), nodes->at(n1),
		nodes->at(n2), nodes->at(n3), nodes->at(n4), E, nu, t));
}
model->add_elements(elements);

// Create restraints
std::vector<Restraint *> *restraints = new std::vector<Restraint *>();
RestraintNode *r1 = new RestraintNode("R1", nodes->at(0));
RestraintNode *r2 = new RestraintNode("R2", nodes->at(static_cast<unsigned long>(N)));
r1->add_all();
r2->add_all();
restraints->push_back(r1);
restraints->push_back(r2);
model->add_restraints(restraints);

// Add distribuited load
std::vector<Load *> *loads = new std::vector<Load *>();
for (int i = 0; i < elements->size(); i++) {
    Membrane *mem = dynamic_cast<Membrane *>(elements->at
        (static_cast<unsigned long>(i)));
    loads->push_back( new LoadMembraneDistributed("DV100kN V @" +
        std::to_string(i + 1), mem, 4, 3, -100, 0, -100, 1));
}
std::vector<LoadPattern *> *loadpattern = new std::vector<LoadPattern *>();
loadpattern->push_back(new LoadPatternConstant("LOADCONSTANT", loads));
model->add_load_patterns(loadpattern);

// Create analysis
StaticAnalysis *analysis = new StaticAnalysis(model);
analysis->analyze(false);

// Save results to file
model->save_results("bridge.txt");

// Delete data
analysis->clear();
delete elements;
delete loadpattern;
delete loads;
delete restraints;
delete nodes;
delete model;
delete analysis;

The results are the following:

If N number of membranes is greater then GPU analysis will be faster.

Author

Pablo Pizarro R. | 2018-2019

License

This project is licensed under MIT [https://opensource.org/licenses/MIT/]