Homec4science
Diffusion LibMultiScale (features/components) Browse

libmultiscale/src/continuum/1d/mesh_meca1d.cc
b05d5f2aba78features/components

Last Change
Raw File

mesh_meca1d.cc
View Options

/**
* @file mesh_meca1d.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Tue Jan 14 17:18:38 2014
*
* @brief The mesh data structure of MECA1D
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include "mesh_meca1d.hh"
#include "ball.hh"
#include "lm_parser.hh"
/* -------------------------------------------------------------------------- */
#include <fstream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
MeshMeca1D::MeshMeca1D() {
}
/* -------------------------------------------------------------------------- */
void MeshMeca1D::buildMesh(Real element_size,Geometry & g,Real mass_density) {
if (element_size < 0) LM_FATAL("element size not set: use ELEM_SIZE keyword");
Ball & b = dynamic_cast<Ball&>(g);
if (element_size == 0) LM_FATAL("element step size should be set by ELEM_SIZE keyword");
UInt nx = (UInt)(nearbyint((b.rMax() - b.rMin())/element_size));
DUMP("nx = " << nx,DBG_MESSAGE);
if(b.rMin() > 0) {
u=(2*nx+2);p0=(2*nx+2);
v=(2*nx+2);f=(2*nx+2);
m=(2*nx+2);//f_old=(2*nx+2);
alpha=(2*nx+2);T=(2*nx+2);Tdt=(2*nx+2);
n_nodes=(2*nx+2);n_elements=(2*nx);
external_heat_rate=(2*nx+2);
heat_rate=(2*nx+2);
elastic_constant = 0.0;
density = mass_density;
elt = NULL;
elt = new FEMeca1D[2*nx];
for (UInt i=0;i< nx;++i) {
elt[i].n1 = i;
elt[i].n2 = i+1;
}
for (UInt i=0;i< nx;++i) {
elt[i+nx].n1 = i+nx+1;
elt[i+nx].n2 = i+nx+2;
}
Real rmin = b.rMin();
Real rmax = b.rMax();
Real center = b.getCenter(0);
Real step = (rmax - rmin)/nx;
e_size = step;
DUMP("taille element " << step,DBG_MESSAGE);
for (UInt i=0; i< nx + 1; ++i) {
u[i] = 0;
p0[i] = (static_cast<Real>(i)-static_cast<Real>(nx))/static_cast<Real>(nx);
v[i] = 0.0;
m[i] = mass_density*step;
f[i] = 0.0;
alpha[i] = 1.0;
T[i] = 0.0;
Tdt[i] = 0.0;
if(i==0 || i==nx)
m[i] /= 2;
}
for (UInt i=0; i< nx + 1; ++i) {
u[i+nx+1] = 0;
p0[i+nx+1] = static_cast<Real>(i)/static_cast<Real>(nx);
v[i+nx+1] = 0.0;
m[i+nx+1] = mass_density*step;
f[i+nx+1] = 0.0;
alpha[i+nx+1] = 1.0;
T[i+nx+1] = 0.0;
Tdt[i+nx+1] = 0.0;
if(i==0 || i == nx)
m[i+nx+1] /= 2;
}
// Scale the nodal positions
for (UInt p=0; p<nx+1; ++p){
p0[p] = p0[p]*(rmax-rmin) - rmin + center;
//p0[p] = p0[p]*(rmax-rmin) - rmin;
DUMP("creating node at position " << p0[p],DBG_DETAIL);
}
for (UInt p=0; p<nx+1; ++p){
p0[p+nx+1] = p0[p+nx+1]*(rmax-rmin) + rmin + center;
//p0[p+nx+1] = p0[p+nx+1]*(rmax-rmin) + rmin;
DUMP("creating node at position " << p0[p+nx+1],DBG_DETAIL);
}
}
else MeshMeca1DWithoutHole(nx,b,mass_density);
}
/* -------------------------------------------------------------------------- */
void MeshMeca1D::MeshMeca1DWithoutHole(UInt nx,Geometry & g,Real mass_density) {
Ball & b = dynamic_cast<Ball&>(g);
u=(2*nx+1);p0=(2*nx+1);
v=(2*nx+1);f=(2*nx+1);
m=(2*nx+1);//f_old=(2*nx+1);
alpha=(2*nx+1);T=(2*nx+1);Tdt=(2*nx+1);
n_nodes=(2*nx+1);n_elements=(2*nx);
external_heat_rate=(2*nx+1);
heat_rate=(2*nx+1);
elastic_constant=0.0;
density = mass_density;
elt = NULL;
elt = new FEMeca1D[2*nx];
for (UInt i=0;i< nx;++i){
elt[i].n1 = i;
elt[i].n2 = i+1;
}
for (UInt i=0;i< nx;++i){
elt[i+nx].n1 = i+nx;
elt[i+nx].n2 = i+nx+1;
}
Real rmin = b.rMin();
Real rmax = b.rMax();
Real center = b.getCenter(0);
Real step = (rmax - rmin)/nx;
e_size = step;
DUMP("element size " << step,DBG_MESSAGE);
for (UInt i=0; i< nx + 1; ++i) {
u[i] = 0;
p0[i] = (static_cast<Real>(i)-static_cast<Real>(nx))/static_cast<Real>(nx);
v[i] = 0.0;
m[i] = mass_density*step;
f[i] = 0.0;
alpha[i] = 1.0;
T[i] = 0.0;
Tdt[i] = 0.0;
if(i==0)
m[i] /= 2;
}
for (UInt i=0; i< nx; ++i){
u[i+nx+1] = 0;
p0[i+nx+1] = static_cast<Real>(i+1)/static_cast<Real>(nx);
v[i+nx+1] = 0.0;
m[i+nx+1] = mass_density*step;
f[i+nx+1] = 0.0;
alpha[i+nx+1] = 1.0;
T[i+nx+1] = 0.0;
Tdt[i+nx+1] = 0.0;
if(i == nx-1)
m[i+nx+1] /= 2;
}
// Scale the nodal positions
for (UInt p=0; p<nx+1; ++p){
p0[p] = p0[p]*(rmax-rmin) - rmin + center;
DUMP("creating node at position " << p0[p],DBG_DETAIL);
}
for (UInt p=0; p<nx; ++p){
p0[p+nx+1] = p0[p+nx+1]*(rmax-rmin) + rmin + center;
DUMP("creating node at position " << p0[p+nx+1],DBG_DETAIL);
}
}
/* -------------------------------------------------------------------------- */
Real MeshMeca1D::readMesh(const std::string & filename,Real mass_density) {
std::ifstream f;
f.open(filename.c_str());
if (!f.is_open()){
LM_FATAL("could not open config file " << filename);
}
UInt line_count = 0;
std::vector<Real> nd_coords;
while (f.good()){
std::string buffer;
std::getline(f,buffer);
++line_count;
unsigned long found = buffer.find_first_of("#");
std::string buffer_comment_free = buffer.substr(0,found-1);
buffer = buffer_comment_free;
// if it is a comment, then do nothing
if (buffer[0] == '#' || buffer == ""){
continue;
}
std::stringstream line(buffer);
Real x;
try {
Parser::parse(x,line);
}
catch (LibMultiScaleException & e) {
LM_THROW(filename << ":" << line_count << ":" <<
" unable to read coordinate: " << buffer);
}
nd_coords.push_back(x);
}
DUMP("Number of read nodes " << nd_coords.size(),DBG_MESSAGE);
UInt n_nodes = nd_coords.size();
this->u = n_nodes;
this->p0 = n_nodes;
this->v = n_nodes;
this->f = n_nodes;
this->m = n_nodes;
// this->f_old = n_nodes;
this->alpha = n_nodes;
this->T = n_nodes;
this->Tdt = n_nodes;
this->n_nodes = n_nodes;
this->external_heat_rate = n_nodes;
this->heat_rate = n_nodes;
this->n_elements = n_nodes-1;
elt = new FEMeca1D[n_elements];
//build connectivity
for (UInt i = 0 ; i < n_elements ; ++i) {
elt[i].n1 = i;
elt[i].n2 = i+1;
}
//build node coords and initial nodal values
for (UInt i = 0; i < n_nodes; ++i) {
this->u[i] = 0;
this->p0[i] = nd_coords[i];
this->v[i] = 0.0;
this->f[i] = 0.0;
this->m[i] = 0.0;
this->alpha[i] = 1.0;
this->T[i] = 0.0;
this->Tdt[i] = 0.0;
}
Real h = -1.;
//assemble manually the mass
for (UInt i = 0; i < n_elements; ++i) {
UInt n1 = elt[i].n1;
UInt n2 = elt[i].n2;
if (h < 0) h = p0[n2] - p0[n1];
if (h != p0[n2] - p0[n1]) LM_FATAL("mesh must be uniform in current state");
m[n1] += mass_density*h/2;
m[n2] += mass_density*h/2;
}
return h;
}
/* -------------------------------------------------------------------------- */
UInt MeshMeca1D::nbElements(){
return n_elements;
}
/* -------------------------------------------------------------------------- */
UInt MeshMeca1D::nbNodes(){
return n_nodes;
}
/* -------------------------------------------------------------------------- */
Real MeshMeca1D::elementSize(){
return e_size;
}
/* -------------------------------------------------------------------------- */
void MeshMeca1D::setE(Real E){
elastic_constant = E;
}
/* -------------------------------------------------------------------------- */
Real MeshMeca1D::getE(){
return elastic_constant;
}
/* -------------------------------------------------------------------------- */
Real MeshMeca1D::getRho(){
return density;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
c4science · Help