1 #ifndef MOAB_PARAMETRIZER_HPP
2 #define MOAB_PARAMETRIZER_HPP
3 #include "moab/Matrix3.hpp"
4 #include "moab/CartVect.hpp"
5 #include "moab/ElemUtil.hpp"
6 namespace moab {
7
8 namespace element_utility {
9 //non-exported functionality
10 namespace {
11
12 template< typename Moab,
13 typename Entity_handle, typename Points>
get_moab_points(Moab & moab,Entity_handle eh,Points & points)14 void get_moab_points( Moab & moab,
15 Entity_handle eh,
16 Points & points){
17 const Entity_handle* connectivity_begin;
18 int num_vertices;
19 moab.get_connectivity( eh, connectivity_begin, num_vertices);
20 //TODO: This is hacky, it works correctly since
21 //CartVect is only double d[ 3], with a default
22 //constructor.. get_coords() should be
23 //flexible enough to allow other types..
24 points.resize( num_vertices);
25 moab.get_coords( connectivity_begin, num_vertices, &(points[ 0][ 0]));
26 }
27
28 } // non-exported functionality
29
30 template< typename Element_map>
31 class Element_parametrizer{
32 public:
33 //public types
34 typedef Element_map Map;
35 private:
36 typedef Element_parametrizer< Map> Self;
37 public: //public functionality
Element_parametrizer()38 Element_parametrizer(): map(){}
Element_parametrizer(const Self & f)39 Element_parametrizer( const Self & f): map( f.map) {}
40 public:
41 template< typename Moab, typename Entity_handle, typename Point>
operator ()(Moab & moab,const Entity_handle & eh,const Point & point,const double tol)42 std::pair< bool, Point> operator()( Moab & moab,
43 const Entity_handle & eh,
44 const Point & point,
45 const double tol){
46 typedef std::vector< moab::CartVect> Points;
47 Points points;
48 get_moab_points( moab, eh, points);
49 return map( moab, eh, points, point, tol);
50 }
51 private:
52 Element_map map;
53 }; //class Element_parametrizer
54
55 class Parametrizer{
56 private:
57 typedef Parametrizer Self;
58 typedef moab::EntityHandle Entity_handle;
59 public: //public functionality
Parametrizer()60 Parametrizer(): hex_map(), tet_map(){}
Parametrizer(const Self & f)61 Parametrizer( const Self & f): hex_map( f.hex_map),
62 tet_map( f.tet_map) {}
63 public:
64 template< typename Moab, typename Entity_handle, typename Point>
operator ()(Moab & moab,const Entity_handle & eh,const Point & point)65 std::pair< bool, Point> operator()( Moab & moab,
66 const Entity_handle & eh,
67 const Point & point){
68 //get entity
69 typedef std::vector< moab::CartVect> Points;
70 Points points;
71 get_moab_points( moab, eh, points);
72 //get type
73 switch( moab.type_from_handle( eh)){
74 case moab::MBHEX:
75 return hex_map( moab, eh,
76 points, point);
77 case moab::MBTET:
78 return tet_map( moab, eh,
79 points, point);
80 //TODO: not correct..
81 //TODO: add quadratic hex, and a proper case for
82 // spectral hex
83 default:
84 quadratic_hex_map( moab, eh,
85 points, point);
86 return spectral_hex_map( moab, eh,
87 points, point);
88 std::cerr << "Element type not supported"
89 << std::endl;
90 return make_pair( false, Point(3, 0.0));
91 }
92 }
93 template< typename Moab, typename Entity_handle, typename Point>
interpolate(Moab & moab,const Entity_handle & eh,const Point & natural_coords)94 void interpolate( Moab & moab, const Entity_handle & eh,
95 const Point & natural_coords){
96 //get field values from moab tag,
97 //then
98 //evaluate_scalar_field();
99 }
100 private:
101 Linear_hex_map< moab::Matrix3> hex_map;
102 Linear_tet_map< Entity_handle, moab::Matrix3> tet_map;
103 Spectral_hex_map< moab::Matrix3> spectral_hex_map;
104 Quadratic_hex_map< moab::Matrix3> quadratic_hex_map;
105 }; //class Parametrizer
106
107 }// namespace element_utility
108 } // namespace moab
109 #endif //MOAB_PARAMETRIZER_HPP
110