Repository navigation
Expand file tree
/
Copy pathHex8_box_utils.hpp
More file actions
174 lines (148 loc) · 5.64 KB
/
Copy pathHex8_box_utils.hpp
File metadata and controls
174 lines (148 loc) · 5.64 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
#ifndef _Hex8_box_utils_hpp_
#define _Hex8_box_utils_hpp_
//@HEADER
// ************************************************************************
//
// HPCCG: Simple Conjugate Gradient Benchmark Code
// Copyright (2006) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER
#include <stdexcept>
#include <box_utils.hpp>
#include <ElemData.hpp>
#include <simple_mesh_description.hpp>
#include <Hex8.hpp>
namespace miniFE {
template<typename GlobalOrdinal>
void get_hex8_node_ids(int nx, int ny,
GlobalOrdinal node0,
GlobalOrdinal* elem_node_ids)
{
//Given box dimensions nx and ny, and a starting node
//(local-node-0 for a hex8), compute the other nodes
//of the hex8 using the exodus ordering convention.
elem_node_ids[0] = node0;
elem_node_ids[1] = node0 + 1;
elem_node_ids[2] = node0 + nx + 1;
elem_node_ids[3] = node0 + nx;
elem_node_ids[4] = node0 + nx*ny;
elem_node_ids[5] = node0 + 1 + nx*ny;
elem_node_ids[6] = node0 + nx + nx*ny + 1;
elem_node_ids[7] = node0 + nx + nx*ny;
}
template<typename Scalar>
void get_hex8_node_coords_3d(Scalar x, Scalar y, Scalar z,
Scalar hx, Scalar hy, Scalar hz,
Scalar* elem_node_coords)
{
//Input: x,y,z are the coordinates of local-node 0 for a Hex8.
//'hx', 'hy', 'hz' are the lengths of the sides of the element
//in each direction.
elem_node_coords[0] = x;
elem_node_coords[1] = y;
elem_node_coords[2] = z;
elem_node_coords[3] = x + hx;
elem_node_coords[4] = y;
elem_node_coords[5] = z;
elem_node_coords[6] = x + hx;
elem_node_coords[7] = y + hy;
elem_node_coords[8] = z;
elem_node_coords[9] = x;
elem_node_coords[10] = y + hy;
elem_node_coords[11] = z;
elem_node_coords[12] = x;
elem_node_coords[13] = y;
elem_node_coords[14] = z + hz;
elem_node_coords[15] = x + hx;
elem_node_coords[16] = y;
elem_node_coords[17] = z + hz;
elem_node_coords[18] = x + hx;
elem_node_coords[19] = y + hy;
elem_node_coords[20] = z + hz;
elem_node_coords[21] = x;
elem_node_coords[22] = y + hy;
elem_node_coords[23] = z + hz;
}
template<typename GlobalOrdinal, typename Scalar>
void
get_elem_nodes_and_coords(const simple_mesh_description<GlobalOrdinal>& mesh,
GlobalOrdinal elemID,
GlobalOrdinal* node_ords, Scalar* node_coords)
{
int global_nodes_x = mesh.global_box[0][1]+1;
int global_nodes_y = mesh.global_box[1][1]+1;
int global_nodes_z = mesh.global_box[2][1]+1;
if (elemID < 0) {
//I don't think this can happen, but check for the sake of paranoia...
throw std::runtime_error("get_elem_nodes_and_coords ERROR, negative elemID");
}
int elem_int_x, elem_int_y, elem_int_z;
get_int_coords(elemID, global_nodes_x-1, global_nodes_y-1, global_nodes_z-1,
elem_int_x, elem_int_y, elem_int_z);
GlobalOrdinal nodeID = get_id<GlobalOrdinal>(global_nodes_x, global_nodes_y, global_nodes_z, elem_int_x, elem_int_y, elem_int_z);
#ifdef MINIFE_DEBUG_VERBOSE
std::cout<<"\nelemID: "<<elemID<<", nodeID: "<<nodeID<<std::endl;
#endif
get_hex8_node_ids(global_nodes_x, global_nodes_y, nodeID, node_ords);
//Map node-IDs to rows because each processor may have a non-contiguous block of
//node-ids, but needs a contiguous block of row-numbers:
#ifdef MINIFE_DEBUG_VERBOSE
std::cout<<"elem "<<elemID<<" nodes: ";
#endif
for(int i=0; i<Hex8::numNodesPerElem; ++i) {
#ifdef MINIFE_DEBUG_VERBOSE
std::cout<<node_ords[i]<<" ";
#endif
node_ords[i] = mesh.map_id_to_row(node_ords[i]);
}
#ifdef MINIFE_DEBUG_VERBOSE
std::cout << std::endl;
#endif
int global_elems_x = mesh.global_box[0][1];
int global_elems_y = mesh.global_box[1][1];
int global_elems_z = mesh.global_box[2][1];
Scalar ix,iy,iz;
get_coords<GlobalOrdinal,Scalar>(nodeID, global_nodes_x,global_nodes_y,global_nodes_z,
ix,iy,iz);
Scalar hx = 1.0/global_elems_x;
Scalar hy = 1.0/global_elems_y;
Scalar hz = 1.0/global_elems_z;
get_hex8_node_coords_3d(ix, iy, iz, hx, hy, hz, node_coords);
#ifdef MINIFE_DEBUG_VERBOSE
int offset = 0;
for(int i=0; i<Hex8::numNodesPerElem; ++i) {
std::cout << "("<<node_coords[offset++]<<","<<node_coords[offset++]<<","<<node_coords[offset++]<<")";
}
std::cout << std::endl;
#endif
}
template<typename GlobalOrdinal, typename Scalar>
void
get_elem_nodes_and_coords(const simple_mesh_description<GlobalOrdinal>& mesh,
GlobalOrdinal elemID,
ElemData<GlobalOrdinal,Scalar>& elem_data)
{
get_elem_nodes_and_coords(mesh, elemID, elem_data.elem_node_ids, elem_data.elem_node_coords);
}
}//namespace miniFE
#endif