Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81102748
diagonal.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Sep 4, 05:27
Size
3 KB
Mime Type
text/x-c
Expires
Fri, Sep 6, 05:27 (2 d)
Engine
blob
Format
Raw Data
Handle
20514745
Attached To
rGOOSEFEM GooseFEM
diagonal.cpp
View Options
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#include <iostream>
#include <Eigen/Eigen>
#include <GooseFEM/GooseFEM.h>
// alias relevant types for <Eigen/Eigen>
using MatD = GooseFEM::MatD;
using MatS = GooseFEM::MatS;
using ColD = GooseFEM::ColD;
using ColS = GooseFEM::ColS;
// =================================================================================================
int main()
{
// --------
// geometry
// --------
// create a mesh
GooseFEM::Mesh::Tri3::Regular mesh ( 2 , 2 );
// element type for extrapolation and quadrature
GooseFEM::Tri3 el;
// extract relevant mesh data
size_t nnode = mesh.nnode(); // number of nodes
size_t nelem = mesh.nelem(); // number of elements
size_t ndim = mesh.ndim (); // number of dimensions (== 2)
size_t nne = mesh.nne (); // number of nodes per element (== 3)
MatS conn = mesh.conn (); // connectivity (nodes of each element) : one row per element
MatD x0 = mesh.coor (); // nodal position in x- and y-direction : one row per node
size_t ndof = nnode * ndim; // total number of DOFs
// DOF-numbers for each vector component of each node (sequential)
MatS dofs = GooseFEM::Mesh::dofs(nnode,ndim);
// perturb position for testing purpose
x0(0,0) -= .5 ; x0(0,1) -= .2 ;
x0(1,0) -= .1 ; x0(1,1) -= .1 ;
x0(2,0) += .2 ; x0(2,1) -= .15;
x0(3,0) += .1 ; x0(3,1) += .2 ;
x0(4,0) += .1 ; x0(4,1) += .1 ;
x0(5,0) -= .2 ; x0(5,1) += .15;
x0(6,0) -= .5 ; x0(6,1) += .4 ;
x0(7,0) -= .1 ; x0(7,1) += .2 ;
x0(8,0) += .2 ; x0(8,1) += .3 ;
// -------------------------------------
// element & quadrature-point - allocate
// -------------------------------------
// element arrays
ColS dof_e(nne*ndim ); // array with DOF numbers of each row/column of "m_e"
MatD x0_e (nne ,ndim); // nodal positions and displacements (column of vectors)
// quadrature-point position and weight factor (this depends on the element type!)
ColD xi_k(nne);
double w = el.QuadGaussWeight()/static_cast<double>(nne);
// ------------------------
// global system - allocate
// ------------------------
// mass matrix: diagonal only
ColD M(ndof);
// density (constant in space)
double rho = 1.0;
// ------------------------
// global system - assemble
// ------------------------
// zero-initialize global system
M.setZero();
// loop over all elements
for ( size_t e = 0 ; e < nelem ; ++e )
{
// - DOF numbers and nodal positions for element "e"
for ( size_t m = 0 ; m < nne ; ++m ) x0_e .row(m) = x0 .row(conn(e,m));
for ( size_t m = 0 ; m < nne ; ++m ) dof_e.segment(m*ndim,ndim) = dofs.row(conn(e,m));
// - loop to integrate and assemble (N.B. quadrature-points coincide with the nodes)
for ( size_t m = 0 ; m < nne ; ++m )
{
// -- define position of the quadrature-point in area coordinates (coincides with node "m")
xi_k.setZero();
xi_k(m) = 1.0;
// -- evaluate the (gradient of the) shape functions
el.eval( x0_e, xi_k, w );
// -- assemble element mass matrix to global system
for ( size_t i = 0; i < ndim ; ++i )
M( dof_e(m*ndim+i) ) += rho * el.V; // N.B. "el.N(m) == 1", so omitted here
}
}
// ---------------
// print to screen
// ---------------
std::cout << M << std::endl;
return 0;
}
Event Timeline
Log In to Comment