Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/dfpm/meshes/axisymmetric_mesh1d.hpp
998871ce7287reactivetransport

Last Change
Raw File

axisymmetric_mesh1d.hpp
View Options

#ifndef SPECMICP_DFPM_MESHES_AXISYMMETRICMESH1D_HPP
#define SPECMICP_DFPM_MESHES_AXISYMMETRICMESH1D_HPP
#include "mesh1d.hpp"
namespace specmicp {
namespace mesh {
//! \brief A generic Axisymmetric 1D mesh
class AxisymmetricMesh1D: public Mesh1D
{
public:
const index_t radius = 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 radius, scalar_t height):
Mesh1D(radius.rows()),
m_height(height),
m_data(radius.rows(), 3),
m_face_radius(radius.rows()-1)
{
// radius
m_data.col(0) = radius;
// Radius of the faces
for (index_t element=0; element<nb_elements(); ++element)
{
m_face_radius(element) = (radius(element+1)+radius(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_PI*height*(std::pow(m_face_radius(node-1), 2)
- std::pow(radius(node), 2));
}
// 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_PI*height*(std::pow(radius(node), 2)
- std::pow(m_face_radius(node), 2));
}
m_data(nb_nodes()-1, cell_vol_1) = 0.0;
}
//! \brief Return the radius of an element
scalar_t get_radius(index_t element, index_t enode) {
return m_data(get_node(element, enode), radius);
}
//! \brief Return the radius of the face
scalar_t get_radius_face(index_t element) {
return m_face_radius(element);
}
//! \brief Return the position (meaning may vary with the mesh) of the node
scalar_t get_position(index_t node) {return m_data(node, radius);}
//! \brief Return the length of an element
scalar_t get_dx(index_t element) {
return get_radius(element, 0) - get_radius(element, 1);}
//! \brief Return the area of the face at the middle of an element
scalar_t get_face_area(index_t element) {
return 2*M_PI*m_face_radius(element)*m_height;}
//! \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_height;
Matrix m_data;
Vector m_face_radius;
};
//! \brief Return a shared_ptr to an axisymmetric mesh
//!
//! \param radius is a vector of the position of the nodes
//! \param height is the height of the sample (or depth)
inline Mesh1DPtr axisymmetric_mesh1d(Vector radius, scalar_t height)
{
return std::static_pointer_cast<Mesh1D>(
std::make_shared<AxisymmetricMesh1D>(radius, height));
}
//! \brief Factory method to build a pointer to a uniform axisymmetric 1D mesh
//!
//! \param nb_nodes number of nodes in the mesh
//! \param radius of the first node (biggest radius in the mesh)
//! \param dx length of an element
//! \param height is the height of the sample (or depth)
//!
//! Note: this method buil an AxisymmetricMesh1D instance while
//! the specmicp::mesh::axisymmetric_uniform_mesh1d method build a
//! specmicp::mesh::AxisymmetricUniformMesh1D instance.
inline Mesh1DPtr uniform_axisymmetric_mesh1d(index_t nb_nodes, scalar_t radius, scalar_t dx, scalar_t height)
{
Vector radius_s(nb_nodes);
for (index_t node=0; node<nb_nodes; ++node)
{
radius_s(node) = radius - node*dx;
}
return std::static_pointer_cast<Mesh1D>(
std::make_shared<AxisymmetricMesh1D>(radius_s, height));
}
} // end namespace mesh
} // end namespace specmicp
#endif // SPECMICP_DFPM_MESHES_AXISYMMETRICMESH1D_HPP
c4science · Help