Branch data Line data Source code
1 : : #include "meshkit/ModelEnt.hpp"
2 : : #include "meshkit/MKCore.hpp"
3 : : #include "moab/Interface.hpp"
4 : : #include "meshkit/FreeSmoothDomain.hpp"
5 : : #include "moab/mesquite/MsqError.hpp"
6 : :
7 : : namespace MeshKit {
8 : :
9 : : using namespace moab;
10 : : using namespace MESQUITE_NS;
11 : :
12 : 2 : FreeSmoothDomain::FreeSmoothDomain( MKCore* core, const MEntVector& entities )
13 : : : MsqCommonIGeom(core->igeom_instance()->instance()),
14 : : haveEntGeomRelTag(false),
15 [ + - ][ + - ]: 2 : moabIface(core->moab_instance())
16 : : {
17 : : // Group input entities and all child entities by dimension.
18 : : // This way if a vertex is contained in all entities in which it
19 : : // is used rather than just the one that owns it, it will still
20 : : // get assigned correctly as long as we work from highest to lowest
21 : : // dimension.
22 [ + - ][ + + ]: 20 : MEntVector ents_by_dim[4], tmp_vec;
[ + - ][ # # ]
23 : 2 : MEntVector::const_iterator i;
24 [ + - ][ + - ]: 4 : for (i = entities.begin(); i != entities.end(); ++i) {
[ + + ]
25 [ + - ]: 2 : ModelEnt* ent = *i;
26 [ + - ]: 2 : int dim = ent->dimension();
27 [ + - ][ - + ]: 2 : if (dim < 0 || dim > 3)
28 [ # # ]: 0 : throw Error(MK_WRONG_DIMENSION, "Entity of invalid dimension %d", dim);
29 [ + - ]: 2 : ents_by_dim[dim].push_back(ent);
30 [ + + ]: 6 : for (int d = 0; d < dim; ++d) {
31 : 4 : tmp_vec.clear();
32 [ + - ]: 4 : ent->get_adjacencies( d, tmp_vec );
33 [ + - ]: 4 : ents_by_dim[d].insert( ents_by_dim[d].end(), tmp_vec.begin(), tmp_vec.end() );
34 : : }
35 : : }
36 : :
37 : : // Create a tag to store geometry handles on mesh entities
38 : : assert( sizeof(iBase_EntityHandle) == sizeof(moab::EntityHandle) );
39 : 2 : const moab::EntityHandle zero = 0;
40 : : moab::ErrorCode rval = moabIface->tag_get_handle( 0,
41 : : 1,
42 : : MB_TYPE_HANDLE,
43 : : entGeomRel,
44 : : moab::MB_TAG_DENSE|moab::MB_TAG_CREAT,
45 [ + - ]: 2 : &zero );
46 [ - + ][ # # ]: 2 : MBERRCHK( rval, moabIface );
[ # # ][ # # ]
47 : 2 : haveEntGeomRelTag = true;
48 : :
49 : : // NOTE: We only go up to a dimension of 2 here!
50 : : // Mesquite only cares about surface elements constrained
51 : : // to a surface and vertices constrained to lie on surfaces,
52 : : // curves, or points.
53 [ + + ]: 8 : for (int d = 0; d < 3; ++d) {
54 [ + - ]: 6 : std::sort( ents_by_dim[d].begin(), ents_by_dim[d].end() );
55 [ + - ]: 6 : MEntVector::iterator end = std::unique( ents_by_dim[d].begin(), ents_by_dim[d].end() );
56 [ + - ][ + - ]: 24 : for (i = ents_by_dim[d].begin(); i != end; ++i) {
[ + - ][ + + ]
57 [ + - ][ + - ]: 18 : EntityHandle set = (*i)->mesh_handle();
58 [ + - ][ + - ]: 18 : iBase_EntityHandle geom = (*i)->geom_handle();
59 [ - + ]: 18 : if (!geom)
60 [ # # ]: 0 : throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Cannot do free smooth if ModelEnt doesn't have geometry" );
61 : :
62 : : // assign surface elements
63 [ + + ]: 18 : if (d == 2) {
64 [ + - ]: 2 : Range elems;
65 [ + - ]: 2 : rval = moabIface->get_entities_by_dimension( set, 2, elems );
66 [ - + ][ # # ]: 2 : MBERRCHK( rval, moabIface );
[ # # ][ # # ]
67 [ + - ]: 2 : rval = moabIface->tag_clear_data( entGeomRel, elems, &geom );
68 [ - + ][ # # ]: 2 : MBERRCHK( rval, moabIface );
[ # # ][ # # ]
69 : : }
70 : :
71 : : // assign vertices
72 [ + - ]: 18 : Range verts;
73 [ + - ]: 18 : rval = moabIface->get_entities_by_dimension( set, 0, verts );
74 [ - + ][ # # ]: 18 : MBERRCHK( rval, moabIface );
[ # # ][ # # ]
75 [ + - ]: 18 : rval = moabIface->tag_clear_data( entGeomRel, verts, &geom );
76 [ - + ][ # # ]: 18 : MBERRCHK( rval, moabIface );
[ # # ][ # # ]
77 : 18 : }
78 [ + + # # ]: 10 : }
79 [ # # ]: 2 : }
80 : :
81 : 4 : FreeSmoothDomain::~FreeSmoothDomain()
82 : : {
83 [ + - ]: 2 : if (haveEntGeomRelTag)
84 : 2 : moabIface->tag_delete( entGeomRel );
85 [ - + ]: 2 : }
86 : :
87 : 8512 : iBase_EntityHandle FreeSmoothDomain::get_geometry( Mesh::EntityHandle mesh_ent ) const
88 : : {
89 : : iBase_EntityHandle result;
90 : 8512 : moab::EntityHandle ent = (moab::EntityHandle)mesh_ent;
91 [ + - ]: 8512 : moab::ErrorCode rval = moabIface->tag_get_data( entGeomRel, &ent, 1, &result );
92 [ - + ][ # # ]: 8512 : MBERRCHK( rval, moabIface );
[ # # ][ # # ]
93 : 8512 : return result;
94 : : }
95 : :
96 : 64 : void FreeSmoothDomain::get_geometry( const Mesh::EntityHandle* mesh_ents,
97 : : size_t num_handles,
98 : : iBase_EntityHandle* geom_ents,
99 : : MsqError& err ) const
100 : : {
101 : 64 : const moab::EntityHandle* ents = (const moab::EntityHandle*)mesh_ents;
102 : 64 : moab::ErrorCode rval = moabIface->tag_get_data( entGeomRel, ents, num_handles, geom_ents );
103 [ - + ]: 64 : if (moab::MB_SUCCESS != rval) {
104 : : MSQ_SETERR(err)(MsqError::INTERNAL_ERROR,
105 [ # # ]: 0 : "Error querying MOAB for geom");
106 : : }
107 : 64 : }
108 : :
109 : :
110 : 2724 : void FreeSmoothDomain::snap_to( Mesh::VertexHandle handle,
111 : : Vector3D& coordinate ) const
112 : : {
113 : 2724 : iBase_EntityHandle geom = get_geometry( handle );
114 [ + - ]: 2724 : if (geom) {
115 : 2724 : int ierr = move_to( geom, coordinate );
116 : 2724 : IBERRCHK(ierr, "iGeom evaluation error");
117 : : }
118 : 2724 : }
119 : :
120 : 4878 : void FreeSmoothDomain::vertex_normal_at( Mesh::VertexHandle handle,
121 : : Vector3D& coordinate ) const
122 : : {
123 : 4878 : iBase_EntityHandle geom = get_geometry( handle );
124 [ + - ]: 4878 : if (geom) {
125 : 4878 : int ierr = normal( geom, coordinate );
126 : 4878 : IBERRCHK(ierr, "iGeom evaluation error");
127 : : }
128 : 4878 : }
129 : :
130 : 4790 : void FreeSmoothDomain::element_normal_at( Mesh::ElementHandle handle,
131 : : Vector3D& coordinate ) const
132 : : {
133 : 4790 : FreeSmoothDomain::vertex_normal_at( handle, coordinate );
134 : 4790 : }
135 : :
136 : 6 : void FreeSmoothDomain::vertex_normal_at( const Mesh::VertexHandle* handle,
137 : : Vector3D coordinates[],
138 : : unsigned count,
139 : : MsqError& err ) const
140 : : {
141 : 6 : tmpHandles.resize( count );
142 : 6 : get_geometry( handle, count, &tmpHandles[0], err );
143 [ - + ][ # # ]: 12 : MSQ_ERRRTN(err);
[ - + ]
144 : :
145 : 6 : int ierr = normal( &tmpHandles[0], coordinates, count );
146 [ - + ]: 6 : if (iBase_SUCCESS != ierr)
147 [ # # ]: 0 : MSQ_SETERR(err)(MsqError::INTERNAL_ERROR,"Failure to evaluate geometry. iBase ErrorCode %d\n", ierr);
148 : : }
149 : :
150 : 58 : void FreeSmoothDomain::domain_DoF( const Mesh::VertexHandle* handle_array,
151 : : unsigned short* dof_array,
152 : : size_t count,
153 : : MsqError& err ) const
154 : : {
155 [ + - ]: 58 : tmpHandles.resize( count );
156 [ + - ][ + - ]: 58 : get_geometry( handle_array, count, &tmpHandles[0], err );
157 [ + - ][ - + ]: 116 : MSQ_ERRRTN(err);
[ # # ][ # # ]
[ - + ]
158 : :
159 : : int ierr, type;
160 [ + + ]: 338 : for (size_t i = 0; i < count; ++i) {
161 [ + - ][ - + ]: 280 : if (!tmpHandles[i])
162 : 0 : dof_array[i] = 3; // don't have geom handles for volumes
163 [ + + ][ + - ]: 280 : else if (i > 0 && tmpHandles[i-1] == tmpHandles[i])
[ + - ][ + + ]
[ + + ]
164 : 58 : dof_array[i] = dof_array[i-1];
165 : : else {
166 [ + - ][ + - ]: 222 : iGeom_getEntType( geomIFace, tmpHandles[i], &type, &ierr );
167 : 222 : dof_array[i] = type;
168 : :
169 [ - + ]: 222 : if (iBase_SUCCESS != ierr) {
170 [ # # ][ # # ]: 0 : MSQ_SETERR(err)(MsqError::INTERNAL_ERROR,"Failure to evaluate geometry. iBase ErrorCode %d\n", ierr);
171 : 0 : return;
172 : : }
173 : : }
174 : : }
175 : : }
176 : :
177 : 910 : void FreeSmoothDomain::closest_point( Mesh::VertexHandle handle,
178 : : const Vector3D& position,
179 : : Vector3D& closest,
180 : : Vector3D& normal,
181 : : MsqError& err ) const
182 : : {
183 : 910 : iBase_EntityHandle geom = get_geometry( handle );
184 [ + - ]: 910 : if (geom) {
185 : 910 : int ierr = closest_and_normal( geom, position, closest, normal );
186 [ - + ]: 910 : if (iBase_SUCCESS != ierr) {
187 [ # # ]: 910 : MSQ_SETERR(err)(MsqError::INTERNAL_ERROR,"Failure to evaluate geometry. iBase ErrorCode %d\n", ierr);
188 : : }
189 : : }
190 : 910 : }
191 : :
192 : :
193 [ + - ][ + - ]: 156 : } // namespace MeshKit
|
