Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected],
25 : : [email protected]
26 : :
27 : : ***************************************************************** */
28 : : #include "PlanarDomain.hpp"
29 : : #include "MsqError.hpp"
30 : : #include "MsqVertex.hpp"
31 : : #include "DomainUtil.hpp"
32 : :
33 : : #include <algorithm>
34 : :
35 [ - + ]: 104 : MBMesquite::PlanarDomain::~PlanarDomain() {}
36 : :
37 : 51 : void MBMesquite::PlanarDomain::set_plane( const MBMesquite::Vector3D& normal, const MBMesquite::Vector3D& point )
38 : : {
39 : 51 : mNormal = normal;
40 : 51 : mNormal.normalize();
41 : 51 : mCoeff = -( mNormal % point );
42 : 51 : }
43 : :
44 : 1421335 : void MBMesquite::PlanarDomain::snap_to( MBMesquite::Mesh::VertexHandle /*entity_handle*/, Vector3D& coordinate ) const
45 : : {
46 [ + - ]: 1421335 : coordinate -= mNormal * ( mNormal % coordinate + mCoeff );
47 : 1421335 : }
48 : :
49 : 0 : void MBMesquite::PlanarDomain::vertex_normal_at( MBMesquite::Mesh::VertexHandle /*entity_handle*/,
50 : : MBMesquite::Vector3D& coordinate ) const
51 : : {
52 : 0 : coordinate = mNormal;
53 : 0 : }
54 : :
55 : 267950 : void MBMesquite::PlanarDomain::element_normal_at( MBMesquite::Mesh::ElementHandle /*entity_handle*/,
56 : : MBMesquite::Vector3D& coordinate ) const
57 : : {
58 : 267950 : coordinate = mNormal;
59 : 267950 : }
60 : :
61 : 96692 : void MBMesquite::PlanarDomain::vertex_normal_at( const MBMesquite::Mesh::VertexHandle*, Vector3D coords[],
62 : : unsigned count, MBMesquite::MsqError& ) const
63 : : {
64 [ + + ]: 784991 : for( unsigned i = 0; i < count; ++i )
65 : 688299 : coords[i] = mNormal;
66 : 96692 : }
67 : :
68 : 229509 : void MBMesquite::PlanarDomain::closest_point( MBMesquite::Mesh::VertexHandle, const MBMesquite::Vector3D& position,
69 : : MBMesquite::Vector3D& closest, MBMesquite::Vector3D& normal,
70 : : MBMesquite::MsqError& ) const
71 : : {
72 : 229509 : normal = mNormal;
73 [ + - ][ + - ]: 229509 : closest = position - mNormal * ( mNormal % position + mCoeff );
74 : 229509 : }
75 : :
76 : 96692 : void MBMesquite::PlanarDomain::domain_DoF( const Mesh::VertexHandle*, unsigned short* dof_array, size_t num_vertices,
77 : : MsqError& ) const
78 : : {
79 [ + - ]: 96692 : std::fill( dof_array, dof_array + num_vertices, 2 );
80 : 96692 : }
81 : :
82 : 4 : void MBMesquite::PlanarDomain::flip()
83 : : {
84 [ + - ]: 4 : mNormal = -mNormal;
85 : 4 : mCoeff = -mCoeff;
86 : 4 : }
87 : :
88 : 4 : void MBMesquite::PlanarDomain::fit_vertices( Mesh* mesh, MsqError& err, double epsilon )
89 : : {
90 : : // Our goal here is to consider only the boundary (fixed) vertices
91 : : // when calculating the plane. If for some reason the user wants
92 : : // to snap a not-quite-planar mesh to a plane during optimization,
93 : : // if possible we want to start with the plane containing the fixed
94 : : // vertices, as those won't get snapped. If no vertices are fixed,
95 : : // then treat them all as fixed for the purpose calculating the plane
96 : : // (consider them all.)
97 : :
98 [ + - ][ + - ]: 8 : std::vector< Mesh::VertexHandle > verts, fixed;
[ + - ]
99 [ + - ][ + - ]: 4 : mesh->get_all_vertices( verts, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
100 [ + - ][ + - ]: 4 : DomainUtil::get_fixed_vertices( mesh, arrptr( verts ), verts.size(), fixed, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
101 : :
102 : 4 : bool do_all_verts = true;
103 [ + - ]: 4 : if( fixed.size() > 2 )
104 : : {
105 : 4 : do_all_verts = false;
106 [ + - ][ + - ]: 4 : fit_vertices( mesh, arrptr( fixed ), fixed.size(), err, epsilon );
107 : :
108 : : // if we failed with only the fixed vertices, try again with all of the
109 : : // vertices
110 [ + - ][ - + ]: 4 : if( err )
111 : : {
112 [ # # ]: 0 : err.clear();
113 : 0 : do_all_verts = true;
114 : : }
115 : : }
116 : :
117 [ - + ]: 4 : if( do_all_verts )
118 : : {
119 [ # # ][ # # ]: 0 : fit_vertices( mesh, arrptr( verts ), verts.size(), err, epsilon );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
120 : : }
121 : :
122 : : // now count inverted elements
123 : 4 : size_t inverted_count = 0, total_count = 0;
124 [ + - ][ + - ]: 8 : std::vector< Mesh::ElementHandle > elems;
125 [ + - ][ + - ]: 8 : std::vector< size_t > junk;
126 [ + - ][ + - ]: 8 : std::vector< MsqVertex > coords;
127 [ + - ][ + - ]: 4 : mesh->get_all_elements( elems, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
128 [ + + ]: 204 : for( size_t i = 0; i < elems.size(); ++i )
129 : : {
130 : :
131 : : EntityTopology type;
132 [ + - ][ + - ]: 200 : mesh->elements_get_topologies( &elems[i], &type, 1, err );
133 [ + - ][ - + ]: 200 : if( TopologyInfo::dimension( type ) != 2 ) continue;
134 : :
135 : 200 : verts.clear();
136 [ + - ][ + - ]: 200 : mesh->elements_get_attached_vertices( arrptr( elems ), 1, verts, junk, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
137 [ - + ]: 200 : if( verts.size() < 3 ) continue;
138 : :
139 [ + - ]: 200 : coords.resize( verts.size() );
140 [ + - ][ + - ]: 200 : mesh->vertices_get_coordinates( arrptr( verts ), arrptr( coords ), 3, err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ - + ]
141 [ + - ][ + - ]: 200 : Vector3D n = ( coords[1] - coords[0] ) * ( coords[2] - coords[0] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
142 : 200 : ++total_count;
143 [ + - ][ + - ]: 200 : if( n % mNormal < 0.0 ) ++inverted_count;
144 : : }
145 : :
146 : : // if most elements are inverted, flip normal
147 [ + - ][ + - ]: 8 : if( 2 * inverted_count > total_count ) this->flip();
[ + - ]
148 : : }
149 : :
150 : 4 : void MBMesquite::PlanarDomain::fit_vertices( MBMesquite::Mesh* mesh, const MBMesquite::Mesh::VertexHandle* verts,
151 : : size_t num_verts, MBMesquite::MsqError& err, double epsilon )
152 : : {
153 [ + - ]: 4 : std::vector< MsqVertex > coords( num_verts );
154 [ + - ][ + - ]: 4 : mesh->vertices_get_coordinates( verts, arrptr( coords ), num_verts, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
155 : :
156 [ + - ][ + - ]: 4 : if( epsilon <= 0.0 ) epsilon = DomainUtil::default_tolerance( arrptr( coords ), num_verts );
[ + - ]
157 : :
158 [ + - ][ + + ]: 16 : Vector3D pts[3];
159 [ + - ][ + - ]: 4 : if( !DomainUtil::non_colinear_vertices( arrptr( coords ), num_verts, pts, epsilon ) )
[ - + ]
160 : : {
161 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "All vertices are colinear", MsqError::INVALID_MESH );
162 : 0 : return;
163 : : }
164 : :
165 [ + - ][ + - ]: 4 : this->set_plane( ( pts[1] - pts[0] ) * ( pts[2] - pts[0] ), pts[0] );
[ + - ][ + - ]
[ + - ]
166 [ + - ][ + - ]: 64 : }
|
