Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/src/dfpm/meshes/generic_mesh1d.hpp
d4f8837fe33areactmicp/unsaturated

Last Change
Raw File

generic_mesh1d.hpp
View Options

#ifndef SPECMICP_DFPM_MESHES_GENERICMESH1D_HPP
#define SPECMICP_DFPM_MESHES_GENERICMESH1D_HPP
#include "mesh1d.hpp"
namespace specmicp {
namespace mesh {
//! \brief A generic 1D mesh
class GenericMesh1D: public Mesh1D
{
public:
const index_t coord = 0;
const index_t cell_vol_0 = 1;
const index_t cell_vol_1 = 2;
//! \brief Build an axisymmetric 1D mesh
//!
//! \param radius is a vector of the position of the nodes
//! \param height is the height of the sample (or depth)
AxisymmetricMesh1D(Vector coords, scalar_t section):
Mesh1D(coords.rows()),
m_section(section),
m_data(radius.rows(), 4),
m_face_coord(coords.rows()-1)
{
// radius
m_data.col(coord) = coords;
// Radius of the faces
for (index_t element=0; element<nb_elements(); ++element)
{
m_face_coord(element) = (coords(element+1)+coords(element))/2.0;
}
// Volume at the left of the cell
m_data(0, cell_vol_0) = 0.0;
for (index_t node=1; node<nb_nodes(); ++node)
{
m_data(node, cell_vol_0) = m_section*(coords(node)-m_face_coord(node-1));
}
// Volume of the cell at the right of the node
for (index_t node=0; node<nb_nodes()-1; ++node)
{
m_data(node, cell_vol_1) = m_section*(m_face_coord(node)-coords(node));
}
m_data(nb_nodes()-1, cell_vol_1) = 0.0;
}
scalar_t get_coord(index_t element, index_t enode) {return m_data(get_node(element, enode), coord);}
//! \brief Return the position (meaning may vary with the mesh) of the node
scalar_t get_position(index_t node) {return m_data(node, coord);}
//! \brief Return the length of an element
scalar_t get_dx(index_t element) {
return get_coord(element, 1) - get_coord(element, 0);}
//! \brief Return the area of the face at the middle of an element
scalar_t get_face_area(index_t element) {
return m_section;}
//! \brief Return the volume of an element
scalar_t get_volume_element(index_t element) {
return (m_data(get_node(element, 0), cell_vol_1)
+ m_data(get_node(element, 1), cell_vol_0));}
//! \brief Return the volume of a cell (element of the dual mesh)
scalar_t get_volume_cell(index_t node) {
return (m_data(node, cell_vol_0) + m_data(node, cell_vol_1));}
//! \brief Return the volume of a cell inside an element
virtual scalar_t get_volume_cell_element(index_t element, index_t enode) {
if (enode == 0)
return m_data(get_node(element, enode), cell_vol_1);
else
return m_data(get_node(element, enode), cell_vol_0);
}
private:
scalar_t m_section;
Matrix m_data;
Vector m_face_coord;
};
//! \brief Return a shared_ptr to an axisymmetric mesh
//!
//! \param coords is a vector of the position of the nodes
//! \param section is the section of the sample (or depth)
inline Mesh1DPtr generic_mesh1d(Vector coords, scalar_t section)
{
return std::static_pointer_cast<Mesh1D>(
std::make_shared<GenericMesh1D>(coords, section));
}
} // end namespace mesh
} // end namespace specmicp
#endif // SPECMICP_DFPM_MESHES_GENERICMESH1D_HPP
c4science · Help