Branch data Line data Source code
1 : : #include "meshkit/PostBL.hpp"
2 : : namespace MeshKit
3 : : {
4 : 1 : int PostBL:: compute_normals()
5 : : // ---------------------------------------------------------------------------
6 : : //! Function: Compute the normal direction for boundary layer compution. \n
7 : : //! Input: From included .hpp file \n
8 : : //! Output: none \n
9 : : // ---------------------------------------------------------------------------
10 : : {
11 : : // COMPUTE NORMALS
12 : : // get tag data and print
13 [ + - ][ + - ]: 1 : all_bl.resize(nodes.size(), 0);
14 : 1 : mb->tag_get_data(MNTag, nodes, &all_bl[0]);
15 : 1 : int count = -1;
16 [ + - ][ + - ]: 9 : for (Range::iterator kter = nodes.begin(); kter != nodes.end(); ++kter){
[ + - ][ + - ]
[ + + ]
17 : 8 : ++count;
18 : : // only one normal in case of single material
19 : 8 : double xdisp = 0.0, ydisp = 0.0, zdisp = 0.0;
20 : :
21 : : // check if node belongs to one or more materials
22 [ + - ][ + - ]: 8 : if(all_bl[count] == 0){
23 [ + - ]: 8 : moab::Range adj_for_normal;
24 [ + - ][ + - ]: 8 : MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_BLDim, false, adj_for_normal, Interface::UNION), mb);
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
25 : : //create the coordinates of the innermost node corresponding to this node
26 [ + - ]: 16 : moab::Range adj_for_norm = intersect(quads, adj_for_normal);
27 : : // now compute the average normal direction for this vertex
28 [ + - ][ + - ]: 8 : moab::CartVect rt(0.0, 0.0, 0.0), v(0.0, 0.0, 0.0);
29 : :
30 [ + - ][ + - ]: 20 : for(Range::iterator qter = adj_for_norm.begin(); qter != adj_for_norm.end(); ++qter){
[ + - ][ + - ]
[ + + ]
31 [ + - ][ + - ]: 12 : MBERRCHK(mb->get_connectivity(&(*qter), 1, adj_qconn),mb);
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
32 : :
33 : 12 : int side_number = 0, sense = 1, offset = 0;
34 [ + - ][ + - ]: 12 : mb->side_number(old_hex[0], (*qter), side_number, sense, offset);
[ + - ]
35 : :
36 [ + - ]: 12 : if(m_GD==3){
37 [ + - ][ + - ]: 12 : get_normal_quad (adj_qconn, v);
38 [ - + ]: 12 : if(sense == 1){
39 : : // do nothing
40 : : }
41 : : else{
42 [ # # ]: 12 : v=-v;
43 : : }
44 : : }
45 [ # # ]: 0 : else if(m_GD==2){
46 [ # # ]: 0 : if(sense == 1){
47 [ # # ][ # # ]: 0 : get_normal_edge(adj_qconn, surf_normal, v);
48 : : }
49 : : else{
50 [ # # ][ # # ]: 0 : get_normal_edge(adj_qconn, -surf_normal, v);
[ # # ]
51 : : }
52 : : }
53 [ + - ]: 12 : rt = rt + v;
54 : : }
55 [ + - ][ + - ]: 8 : if(rt.length() !=0){
56 [ + - ][ + - ]: 8 : xdisp=rt[0]/rt.length();
57 [ + - ][ + - ]: 8 : ydisp=rt[1]/rt.length();
58 [ + - ][ + - ]: 8 : zdisp=rt[2]/rt.length();
59 : : }
60 : : else{
61 : 0 : xdisp=0.0;
62 : 0 : ydisp=0.0;
63 : 0 : zdisp=0.0;
64 : 8 : }
65 : : }
66 [ # # ][ # # ]: 0 : else if(all_bl[count] > 0 && fixmat != -1){ // node belongs to more than one material and fixmat specified
[ # # ][ # # ]
67 : : // NODE B/W MATERIALS
68 [ # # ]: 0 : moab::Range adj_for_normal;
69 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_adjacencies(&(*kter), 1, m_BLDim, false, adj_for_normal, Interface::UNION), mb);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
70 : : //create the coordinates of the innermost node corresponding to this node
71 [ # # ]: 0 : moab::Range adj_for_norm = intersect(quads, adj_for_normal);
72 : : // now compute the average normal direction for this vertex
73 [ # # ][ # # ]: 0 : moab::CartVect rt(0.0, 0.0, 0.0), v(0.0, 0.0, 0.0);
74 : :
75 [ # # ][ # # ]: 0 : for(Range::iterator qter = adj_for_norm.begin(); qter != adj_for_norm.end(); ++qter){
[ # # ][ # # ]
[ # # ]
76 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_connectivity(&(*qter), 1, adj_qconn),mb);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
77 : :
78 [ # # ]: 0 : moab::Range this_quad_hex;
79 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_adjacencies(&(*qter), 1, m_GD, false, this_quad_hex, moab::Interface::UNION),mb);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
80 [ # # ]: 0 : moab::Range quad_hex = intersect(fixmat_ents, this_quad_hex);
81 : :
82 : 0 : int side_number = 0, sense = 1, offset = 0;
83 [ # # ]: 0 : Range::iterator hexter = quad_hex.begin();
84 [ # # ][ # # ]: 0 : mb->side_number(*hexter, (*qter), side_number, sense, offset);
[ # # ]
85 : :
86 [ # # ]: 0 : if(m_GD==3){
87 [ # # ][ # # ]: 0 : get_normal_quad (adj_qconn, v);
88 [ # # ]: 0 : if(sense == 1 ){
89 [ # # ]: 0 : v=-v;
90 : : }
91 : : else{
92 : : // do nothing
93 : : }
94 : : }
95 [ # # ]: 0 : else if(m_GD==2){
96 [ # # ]: 0 : if(sense == 1 ){
97 [ # # ][ # # ]: 0 : get_normal_edge(adj_qconn, -surf_normal, v);
[ # # ]
98 : : }
99 : : else{
100 [ # # ][ # # ]: 0 : get_normal_edge(adj_qconn, surf_normal, v);
101 : : }
102 : : }
103 [ # # ]: 0 : rt = rt + v;
104 : 0 : }
105 [ # # ][ # # ]: 0 : if(rt.length() !=0){
106 [ # # ][ # # ]: 0 : xdisp=rt[0]/rt.length();
107 [ # # ][ # # ]: 0 : ydisp=rt[1]/rt.length();
108 [ # # ][ # # ]: 0 : zdisp=rt[2]/rt.length();
109 : : }
110 : : else{
111 : 0 : xdisp=0.0;
112 : 0 : ydisp=0.0;
113 : 0 : zdisp=0.0;
114 : 0 : }
115 : : }
116 [ # # ][ # # ]: 0 : else if(all_bl[count] > 0 && fixmat == -1){ // node belongs to more than one material and fixmat not specified
[ # # ][ # # ]
117 : : // NODE ON BOUNDARY
118 : : // get the edge that is not on the boundary and count how many such edges we have
119 [ # # ]: 0 : moab::Range adj_for_normal;
120 : 0 : int nEdgeDim = 1;
121 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_adjacencies(&(*kter), 1, nEdgeDim, true, adj_for_normal, Interface::UNION), mb);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
122 [ # # ]: 0 : moab::Range edge_normal;
123 [ # # ]: 0 : if(m_GD == 2)
124 [ # # ][ # # ]: 0 : edge_normal = subtract(adj_for_normal, quads);
125 [ # # ]: 0 : else if(m_GD == 3)
126 [ # # ][ # # ]: 0 : edge_normal = subtract(adj_for_normal, edges);
127 : :
128 : :
129 [ # # ][ # # ]: 0 : if(edge_normal.size() > 1){
130 : : double ncoord[3];
131 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_coords(&(*kter), 1, ncoord),mb);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
132 : :
133 [ # # ][ # # ]: 0 : m_LogFile << "Multiple normals are needed at: " << ncoord[0]
134 [ # # ][ # # ]: 0 : << ", " << ncoord[1] << ", " << ncoord[2] << " #normals " << edge_normal.size() << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
135 : 0 : exit(0);
136 : : }
137 : : else{
138 [ # # ][ # # ]: 0 : m_LogFile << "We've one edge seperating materials 1, edge normal is needed" << edge_normal.size() << std::endl;
[ # # ][ # # ]
139 [ # # ]: 0 : moab::Range edge_conn;
140 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_connectivity(&(*edge_normal.begin()), 1, edge_conn),mb);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
141 : : // now get the normal direction for this edge
142 [ # # ][ # # ]: 0 : moab::CartVect coords[2], bl_coord[1];
[ # # ][ # # ]
143 [ # # ]: 0 : moab::CartVect AB;
144 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_coords(&(*kter), 1, (double*) &bl_coord[0]), mb);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
145 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_coords(edge_conn, (double*) &coords[0]), mb);
[ # # ][ # # ]
[ # # ][ # # ]
146 [ # # ]: 0 : for(int d = 0; d<2; d++){
147 [ # # ][ # # ]: 0 : if(bl_coord[0][0] == coords[d][0] &&
[ # # ][ # # ]
148 [ # # ][ # # ]: 0 : bl_coord[0][1] == coords[d][1] &&
[ # # ][ # # ]
149 [ # # ][ # # ]: 0 : bl_coord[0][2] == coords[d][2]){
150 : : // do nothing
151 : : }
152 : : else{
153 [ # # ]: 0 : AB = (bl_coord[0] - coords[d]);
154 : : }
155 : : }
156 : :
157 [ # # ][ # # ]: 0 : xdisp=AB[0]/AB.length();
158 [ # # ][ # # ]: 0 : ydisp=AB[1]/AB.length();
159 [ # # ][ # # ]: 0 : zdisp=AB[2]/AB.length();
160 : 0 : }
161 : : }
162 [ # # ][ # # ]: 0 : else if(all_bl[count] < 0 ){
163 [ # # ][ # # ]: 0 : m_LogFile << "Material must have an associated BLNode: Error! - shouldn't have gotten here: " << count << std::endl;
[ # # ]
164 : : }
165 : :
166 : : // after the normal computation is done create new BL nodes
167 : 8 : double coords_bl_quad[3], move = 0.0;
168 [ + - ][ + - ]: 8 : MBERRCHK(mb->get_coords(&(*kter), 1, coords_bl_quad),mb);
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
169 : :
170 : 8 : double temp = 0;
171 : 8 : double num = m_Thickness*(m_Bias-1)*(pow(m_Bias, m_Intervals -1));
172 : 8 : double deno = pow(m_Bias, m_Intervals) - 1;
173 [ + - ]: 8 : if (deno !=0)
174 : 8 : temp = num/deno;
175 : : else
176 : 0 : temp = m_Thickness/m_Intervals;
177 : :
178 : : // loop thru intervals to create BL nodes
179 [ + + ]: 24 : for(int j=0; j< m_Intervals; j++){
180 : :
181 : 16 : move+= temp/pow(m_Bias,j);
182 [ - + ]: 16 : if (debug){
183 [ # # ][ # # ]: 0 : m_LogFile << " move:" << move;
184 : : }
185 : : // now compute the coords of the new vertex
186 : 16 : coords_new_quad[0] = coords_bl_quad[0]-move*xdisp;
187 : 16 : coords_new_quad[1] = coords_bl_quad[1]-move*ydisp;
188 : 16 : coords_new_quad[2] = coords_bl_quad[2]-move*zdisp;
189 : :
190 : 16 : int nid = count*m_Intervals+j;
191 : : // Possible TODO's
192 : : //TODO: Check if this vertex is possible (detect possible collision with geometry)
193 : : // TODO: See posibility of using ray tracing
194 : : // TODO: Parallize: Profile T-junction model and try to device an algorithm
195 : : // TODO: Modularize node creation part and use doxygen for all code and design of code, python design and test cases - current functions in code:
196 : : // Setup this, Execute this -- break info sub functions and classes,
197 : : // prepareIO --make this optional when using python,
198 : : // get normal (2d and 3d) -- can be combined to one function
199 : : // get det jacobian (hex elements) --needs check for other elements
200 : : //
201 [ + - ][ + - ]: 16 : MBERRCHK(mb->create_vertex(coords_new_quad, new_vert[nid]), mb);
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
202 [ - + ]: 16 : if (debug){
203 [ # # ][ # # ]: 0 : m_LogFile << std::setprecision (3) << std::scientific << " : created node:" << (nid + 1)
[ # # ][ # # ]
[ # # ]
204 [ # # ][ # # ]: 0 : << " of " << new_vert.size() << " new nodes:- coords: " << coords_new_quad[0]
[ # # ][ # # ]
205 [ # # ][ # # ]: 0 : << ", " << coords_new_quad[1] << ", " << coords_new_quad[2] << std::endl;
[ # # ][ # # ]
[ # # ]
206 : : }
207 : : }
208 : : }
209 : :
210 : 1 : return 0;
211 : : }
212 : :
213 : 12 : void PostBL::get_normal_quad (std::vector<EntityHandle>conn, moab::CartVect &v)
214 : : // ---------------------------------------------------------------------------
215 : : //! Function: Get normal of a quad \n
216 : : //! Input: conn \n
217 : : //! Output: moab::CartVect v \n
218 : : // ---------------------------------------------------------------------------
219 : : {
220 [ + - ][ + + ]: 48 : moab::CartVect coords[3];
221 [ + - ][ + - ]: 12 : MBERRCHK(mb->get_coords(&conn[0], 3, (double*) &coords[0]), mb);
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
222 [ + - ]: 12 : moab::CartVect AB(coords[1] - coords[0]);
223 [ + - ]: 12 : moab::CartVect BC(coords[2] - coords[1]);
224 [ + - ]: 12 : moab::CartVect normal = AB*BC;
225 [ + - ][ + - ]: 12 : normal = normal/normal.length();
226 : 12 : v = normal;
227 : 12 : }
228 : :
229 : 0 : void PostBL::get_normal_edge (std::vector<EntityHandle>conn, moab::CartVect BC, moab::CartVect &v)
230 : : // ---------------------------------------------------------------------------
231 : : //! Function: Get normal of a edge along its quad \n
232 : : //! Input: conn of edge, normal to the surf \n
233 : : //! Output: moab::CartVect v \n
234 : : // ---------------------------------------------------------------------------
235 : : {
236 [ # # ][ # # ]: 0 : moab::CartVect coords[2];
237 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_coords(&conn[0], 2, (double*) &coords[0]), mb);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
238 [ # # ]: 0 : moab::CartVect AB(coords[1] - coords[0]);
239 [ # # ]: 0 : moab::CartVect normal = AB*BC;
240 [ # # ][ # # ]: 0 : normal = normal/normal.length();
241 : 0 : v = normal;
242 : 0 : }
243 : :
244 : 0 : void PostBL::get_det_jacobian(std::vector<moab::EntityHandle> conn, int offset, double &AvgJ)
245 : : // ---------------------------------------------------------------------------
246 : : //! Function: Get determinant of jacobian \n
247 : : //! Input: conn \n
248 : : //! Output: vector x, y and z \n
249 : : // ---------------------------------------------------------------------------
250 : : {
251 : : //TODO: Add quality check for tri/quad and pyramids
252 [ # # ]: 0 : if(m_Conn ==8){
253 : 0 : ++m_JacCalls;
254 [ # # ][ # # ]: 0 : moab::CartVect vertex[8], xi;
[ # # ]
255 [ # # ]: 0 : mstream m_LogFile;
256 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_coords(&conn[offset], 8, (double*) &vertex[0]), mb);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
257 : :
258 : : double corner[8][3] = { { -1, -1, -1 },
259 : : { 1, -1, -1 },
260 : : { 1, 1, -1 },
261 : : { -1, 1, -1 },
262 : : { -1, -1, 1 },
263 : : { 1, -1, 1 },
264 : : { 1, 1, 1 },
265 : 0 : { -1, 1, 1 } };
266 : :
267 [ # # ]: 0 : for (unsigned j = 0; j < 8; ++j) {
268 [ # # ]: 0 : xi[0] = corner[j][0];
269 [ # # ]: 0 : xi[1] = corner[j][1];
270 [ # # ]: 0 : xi[2] = corner[j][2];
271 [ # # ]: 0 : Matrix3 J(0.0);
272 : 0 : double detJ = 0;
273 [ # # ]: 0 : for (unsigned i = 0; i < 8; ++i) {
274 [ # # ]: 0 : const double xi_p = 1 + xi[0]*corner[i][0];
275 [ # # ]: 0 : const double eta_p = 1 + xi[1]*corner[i][1];
276 [ # # ]: 0 : const double zeta_p = 1 + xi[2]*corner[i][2];
277 : 0 : const double dNi_dxi = corner[i][0] * eta_p * zeta_p;
278 : 0 : const double dNi_deta = corner[i][1] * xi_p * zeta_p;
279 : 0 : const double dNi_dzeta = corner[i][2] * xi_p * eta_p;
280 [ # # ][ # # ]: 0 : J(0,0) += dNi_dxi * vertex[i][0];
281 [ # # ][ # # ]: 0 : J(1,0) += dNi_dxi * vertex[i][1];
282 [ # # ][ # # ]: 0 : J(2,0) += dNi_dxi * vertex[i][2];
283 [ # # ][ # # ]: 0 : J(0,1) += dNi_deta * vertex[i][0];
284 [ # # ][ # # ]: 0 : J(1,1) += dNi_deta * vertex[i][1];
285 [ # # ][ # # ]: 0 : J(2,1) += dNi_deta * vertex[i][2];
286 [ # # ][ # # ]: 0 : J(0,2) += dNi_dzeta * vertex[i][0];
287 [ # # ][ # # ]: 0 : J(1,2) += dNi_dzeta * vertex[i][1];
288 [ # # ][ # # ]: 0 : J(2,2) += dNi_dzeta * vertex[i][2];
289 : : }
290 [ # # ]: 0 : J *= 0.125;
291 [ # # ]: 0 : detJ = J.determinant();
292 [ # # ]: 0 : if(detJ <= 0.0){
293 [ # # ][ # # ]: 0 : m_LogFile << "We've negative jacobian at the hex corner: "<< j+1 << std::endl;
[ # # ]
294 : 0 : exit(0);
295 : : }
296 : 0 : AvgJ+=detJ;
297 : : }
298 : 0 : AvgJ/=8;
299 [ # # ]: 0 : if(m_JacCalls == 1){
300 : 0 : m_JLo = AvgJ;
301 : 0 : m_JHi = AvgJ;
302 : : }
303 [ # # ]: 0 : else if(AvgJ < m_JLo){
304 : 0 : m_JLo = AvgJ;
305 : : }
306 [ # # ]: 0 : else if(AvgJ > m_JHi){
307 : 0 : m_JHi = AvgJ;
308 : 0 : }
309 : : }
310 : 0 : }
311 : 0 : void PostBL::find_min_edge_length(moab::Range adj_qn, moab::EntityHandle node , moab::Range bl_nodes, double &e_len)
312 : : // ---------------------------------------------------------------------------
313 : : //! Function: Get minimum edge length from several adjacent quads/edges specified \n
314 : : //! Input: verts of quads, BL vert \n
315 : : //! Output: distance b/w BL vert and inner vert \n
316 : : // ---------------------------------------------------------------------------
317 : : {
318 : : // get nodes adj(a) to BL node
319 [ # # ][ # # ]: 0 : moab::CartVect coords[1];
320 : 0 : double len = 0; e_len = 0;
321 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_coords(&node, 1, (double*) &coords[0]), mb);
[ # # ][ # # ]
[ # # ][ # # ]
322 [ # # ]: 0 : moab::Range non_bl;
323 [ # # ][ # # ]: 0 : non_bl = subtract(adj_qn, bl_nodes);
324 [ # # ]: 0 : int n_non_bl = (int) non_bl.size();
325 : :
326 : : // if there are no bl nodes - this case has already been dealt with
327 [ # # ][ # # ]: 0 : if(non_bl.size() > 0){
328 [ # # ][ # # ]: 0 : moab::CartVect non_bl_coords[4];
329 [ # # ]: 0 : for(int i=0; i< n_non_bl; i++){
330 [ # # ][ # # ]: 0 : MBERRCHK(mb->get_coords(non_bl,(double*) &non_bl_coords[0]), mb);
[ # # ][ # # ]
[ # # ][ # # ]
331 [ # # ]: 0 : moab::CartVect edge(coords[0] - non_bl_coords[0]);
332 [ # # ]: 0 : len = edge.length();
333 [ # # ]: 0 : if(i==0)
334 : 0 : e_len = len;
335 [ # # ]: 0 : if(len < e_len)
336 : 0 : e_len = len;
337 : : }
338 : 0 : }
339 : : // m_LogFile << " node minimum edge length" << e_len << std::endl;
340 : 0 : }
341 [ + - ][ + - ]: 156 : }
|
