Page MenuHomec4science
Files F91614329

axisymmetric_uniform_mesh1d.hpp
No OneTemporary
Actions

File Metadata

Created
Tue, Nov 12, 17:49

axisymmetric_uniform_mesh1d.hpp
View Options

/*-------------------------------------------------------
- Module : reactmicp/meshes
- File : axisymmetric_uniform_mesh1d.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_REACTMICP_MESHES_AXISYMMETRICMESH1D_HPP
#define SPECMICP_REACTMICP_MESHES_AXISYMMETRICMESH1D_HPP
#include "mesh1d.hpp"
// \file axisymmetric_uniform_mesh1d.hpp Uniform axisymmetric 1D mesh
namespace specmicp {
namespace mesh {
//! \brief A uniform 1D mesh
class AxisymmetricUniformMesh1D: public Mesh1D
{
public:
AxisymmetricUniformMesh1D(index_t nb_nodes, scalar_t radius, scalar_t height):
Mesh1D(nb_nodes),
m_radius(radius),
m_radius_int(0.0),
m_dx(radius/(nb_nodes-1)),
m_height(height)
{}
AxisymmetricUniformMesh1D(index_t nb_nodes, scalar_t radius, scalar_t dx, scalar_t height):
Mesh1D(nb_nodes),
m_radius(radius),
m_radius_int(radius-(nb_nodes-1)*dx),
m_dx(dx),
m_height(height)
{}
scalar_t get_radius_node(index_t node) {
return m_radius-m_dx*node;
}
//! Return the radius of a face of an element
scalar_t get_radius_face(index_t element) {
return get_radius_node(get_node(element, 1))+m_dx/2;
}
scalar_t get_position(index_t node) {return get_radius_node(node);}
scalar_t get_dx(index_t _) {return m_dx;}
scalar_t get_face_area(index_t element) {
return 2*M_PI*get_radius_face(element)*m_height;}
scalar_t get_volume_element(index_t element) {return M_PI*m_height*(
std::pow(get_radius_node(get_node(element, 0)),2)
- std::pow(get_radius_node(get_node(element, 1)),2))
;}
scalar_t get_volume_cell(index_t node) {
if (node ==0)
return M_PI*m_height*(std::pow(m_radius,2)-std::pow(get_radius_face(0),2));
else if (node == nb_nodes()-1)
if (m_radius_int == 0.0)
return M_PI*m_height*std::pow(m_dx/2.0,2);
else
return M_PI*m_height*(std::pow(get_radius_face(node-1),2) - std::pow(m_radius_int,2));
else
return M_PI*m_height*(std::pow(get_radius_face(node-1),2) - std::pow(get_radius_face(node),2));
}
scalar_t get_volume_cell_element(index_t element, index_t enode) {
if (enode == 0)
{
return M_PI*m_height*(std::pow(get_radius_node(get_node(element,enode)),2)
- std::pow(get_radius_face(element),2));
}
else
{
return M_PI*m_height*(std::pow(get_radius_face(element),2)
- std::pow(get_radius_node(get_node(element,enode)),2));
}
}
private:
scalar_t m_radius;
scalar_t m_radius_int;
scalar_t m_dx;
scalar_t m_height;
};
//! \brief Factory method to build a pointer to a uniform 1D mesh
inline Mesh1DPtr axisymmetric_uniform_mesh1d(index_t nb_nodes, scalar_t radius, scalar_t height)
{
return std::static_pointer_cast<Mesh1D>(
std::make_shared<AxisymmetricUniformMesh1D>(nb_nodes, radius, height));
}
//! \brief Factory method to build a pointer to a uniform 1D mesh
inline Mesh1DPtr axisymmetric_uniform_mesh1d(index_t nb_nodes, scalar_t radius, scalar_t dx, scalar_t height)
{
return std::static_pointer_cast<Mesh1D>(
std::make_shared<AxisymmetricUniformMesh1D>(nb_nodes, radius, dx, height));
}
} // end namespace mesh
} // end namespace specmicp
#endif // SPECMICP_REACTMICP_MESHES_AXISYMMETRICMESH1D_HPP

Event Timeline

c4science · Help