Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F71393082
mesh_meca1d.cc
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
Thu, Jul 11, 11:45
Size
9 KB
Mime Type
text/x-c
Expires
Sat, Jul 13, 11:45 (2 d)
Engine
blob
Format
Raw Data
Handle
18945761
Attached To
rLIBMULTISCALE LibMultiScale
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 "mesh_meca1d.hh"
#include "ball.hh"
#include "lm_common.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.resize(2 * nx + 2, 1);
p0.resize(2 * nx + 2, 1);
v.resize(2 * nx + 2, 1);
f.resize(2 * nx + 2, 1);
a.resize(2 * nx + 2, 1);
m.resize(2 * nx + 2, 1);
// 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.resize(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.resize(2 * nx + 1, 1);
p0.resize(2 * nx + 1, 1);
v.resize(2 * nx + 1, 1);
f.resize(2 * nx + 1, 1);
a.resize(2 * nx + 1, 1);
m.resize(2 * nx + 1, 1);
// alpha.resize(2 * nx + 1, 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.resize(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.resize(n_nodes, 1);
this->p0.resize(n_nodes, 1);
this->v.resize(n_nodes, 1);
this->f.resize(n_nodes, 1);
this->a.resize(n_nodes, 1);
this->m.resize(n_nodes, 1);
// 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.resize(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__
Event Timeline
Log In to Comment