FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
hecmw_couple_judge.c
Go to the documentation of this file.
1 /*****************************************************************************
2  * Copyright (c) 2019 FrontISTR Commons
3  * This software is released under the MIT License, see LICENSE.txt
4  *****************************************************************************/
5 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <string.h>
9 #include <assert.h>
10 #include <errno.h>
11 #include <math.h>
12 
13 #include "hecmw_struct.h"
14 
15 #include "hecmw_couple_define.h"
16 #include "hecmw_couple_struct.h"
17 #include "hecmw_couple_judge.h"
18 
19 #define FRAC_1_2 (0.5)
20 
21 #define FRAC_1_3 (0.33333333333333333)
22 
23 #define FRAC_1_4 (0.25)
24 
25 #define EPS_PRODUCT (1.0E-06)
26 #define EPS_DISTANCE (1.0E-02)
27 
28 /*================================================================================================*/
29 
30 static void calc_gravity_tri(double p1_x, double p1_y, double p1_z, double p2_x,
31  double p2_y, double p2_z, double p3_x, double p3_y,
32  double p3_z, double *g_x, double *g_y,
33  double *g_z) {
34  *g_x = (p1_x + p2_x + p3_x) * FRAC_1_3;
35  *g_y = (p1_y + p2_y + p3_y) * FRAC_1_3;
36  *g_z = (p1_z + p2_z + p3_z) * FRAC_1_3;
37 }
38 
39 static void calc_gravity_quad(double p1_x, double p1_y, double p1_z,
40  double p2_x, double p2_y, double p2_z,
41  double p3_x, double p3_y, double p3_z,
42  double p4_x, double p4_y, double p4_z,
43  double *g_x, double *g_y, double *g_z) {
44  *g_x = (p1_x + p2_x + p3_x + p4_x) * FRAC_1_4;
45  *g_y = (p1_y + p2_y + p3_y + p4_y) * FRAC_1_4;
46  *g_z = (p1_z + p2_z + p3_z + p4_z) * FRAC_1_4;
47 }
48 
49 static void calc_normal_vector_tri(double p1_x, double p1_y, double p1_z,
50  double p2_x, double p2_y, double p2_z,
51  double p3_x, double p3_y, double p3_z,
52  double *vn_x, double *vn_y, double *vn_z) {
53  double v1_x, v1_y, v1_z;
54  double v2_x, v2_y, v2_z;
55 
56  v1_x = p2_x - p1_x;
57  v1_y = p2_y - p1_y;
58  v1_z = p2_z - p1_z;
59 
60  v2_x = p3_x - p1_x;
61  v2_y = p3_y - p1_y;
62  v2_z = p3_z - p1_z;
63 
64  *vn_x = v1_y * v2_z - v1_z * v2_y;
65  *vn_y = v1_z * v2_x - v1_x * v2_z;
66  *vn_z = v1_x * v2_y - v1_y * v2_x;
67 }
68 
69 static void calc_normal_vector_quad(double p1_x, double p1_y, double p1_z,
70  double p2_x, double p2_y, double p2_z,
71  double p3_x, double p3_y, double p3_z,
72  double p4_x, double p4_y, double p4_z,
73  double *vn_x, double *vn_y, double *vn_z) {
74  double vn1_x, vn1_y, vn1_z;
75  double vn2_x, vn2_y, vn2_z;
76 
77  calc_normal_vector_tri(p1_x, p1_y, p1_z, p2_x, p2_y, p2_z, p4_x, p4_y, p4_z,
78  &vn1_x, &vn1_y, &vn1_z);
79  calc_normal_vector_tri(p3_x, p3_y, p3_z, p4_x, p4_y, p4_z, p2_x, p2_y, p2_z,
80  &vn2_x, &vn2_y, &vn2_z);
81 
82  *vn_x = (vn1_x + vn2_x) * FRAC_1_2;
83  *vn_y = (vn1_y + vn2_y) * FRAC_1_2;
84  *vn_z = (vn1_z + vn2_z) * FRAC_1_2;
85 }
86 
87 static double calc_dot_product(double v1_x, double v1_y, double v1_z,
88  double v2_x, double v2_y, double v2_z,
89  double *dot_product) {
90  *dot_product = v1_x * v2_x + v1_y * v2_y + v1_z * v2_z;
91 
92  return *dot_product;
93 }
94 
95 static int calc_dot_product_tri(double p1_x, double p1_y, double p1_z,
96  double p2_x, double p2_y, double p2_z,
97  double p3_x, double p3_y, double p3_z,
98  double p_x, double p_y, double p_z,
99  double *dot_product, double *distance) {
100  double vn_x, vn_y, vn_z;
101  double g_x, g_y, g_z;
102  double vp_x, vp_y, vp_z;
103 
104  calc_normal_vector_tri(p1_x, p1_y, p1_z, p2_x, p2_y, p2_z, p3_x, p3_y, p3_z,
105  &vn_x, &vn_y, &vn_z);
106  calc_gravity_tri(p1_x, p1_y, p1_z, p2_x, p2_y, p2_z, p3_x, p3_y, p3_z, &g_x,
107  &g_y, &g_z);
108 
109  vp_x = p_x - g_x;
110  vp_y = p_y - g_y;
111  vp_z = p_z - g_z;
112 
113  calc_dot_product(vn_x, vn_y, vn_z, vp_x, vp_y, vp_z, dot_product);
114 
115  *distance =
116  fabs(*dot_product / sqrt(vn_x * vn_x + vn_y * vn_y + vn_z * vn_z));
117 
118  return *dot_product;
119 }
120 
121 static int calc_dot_product_quad(double p1_x, double p1_y, double p1_z,
122  double p2_x, double p2_y, double p2_z,
123  double p3_x, double p3_y, double p3_z,
124  double p4_x, double p4_y, double p4_z,
125  double p_x, double p_y, double p_z,
126  double *dot_product, double *distance) {
127  double vn_x, vn_y, vn_z;
128  double g_x, g_y, g_z;
129  double vp_x, vp_y, vp_z;
130 
131  calc_normal_vector_quad(p1_x, p1_y, p1_z, p2_x, p2_y, p2_z, p3_x, p3_y, p3_z,
132  p4_x, p4_y, p4_z, &vn_x, &vn_y, &vn_z);
133 
134  calc_gravity_quad(p1_x, p1_y, p1_z, p2_x, p2_y, p2_z, p3_x, p3_y, p3_z, p4_x,
135  p4_y, p4_z, &g_x, &g_y, &g_z);
136 
137  vp_x = p_x - g_x;
138  vp_y = p_y - g_y;
139  vp_z = p_z - g_z;
140 
141  calc_dot_product(vn_x, vn_y, vn_z, vp_x, vp_y, vp_z, dot_product);
142 
143  *distance =
144  fabs(*dot_product / sqrt(vn_x * vn_x + vn_y * vn_y + vn_z * vn_z));
145 
146  return *dot_product;
147 }
148 
149 extern int HECMW_couple_judge_tet1(const struct hecmwST_local_mesh *local_mesh,
150  int elem, int surf_id, double coord_px,
151  double coord_py, double coord_pz,
152  double *dot_product, double *distance) {
153  double coord_x[4], coord_y[4], coord_z[4], _dot_product[4], _distance[4];
154  int node_index, node_id, n_positive;
155  int i;
156 
157  node_index = local_mesh->elem_node_index[elem - 1];
158  for (i = 0; i < 4; i++) {
159  node_id = local_mesh->elem_node_item[node_index + i];
160  coord_x[i] = local_mesh->node[3 * (node_id - 1)];
161  coord_y[i] = local_mesh->node[3 * (node_id - 1) + 1];
162  coord_z[i] = local_mesh->node[3 * (node_id - 1) + 2];
163  }
164 
165  calc_dot_product_tri(coord_x[1], coord_y[1], coord_z[1], coord_x[2],
166  coord_y[2], coord_z[2], coord_x[3], coord_y[3],
167  coord_z[3], coord_px, coord_py, coord_pz,
168  &_dot_product[0], &_distance[0]);
169  calc_dot_product_tri(coord_x[0], coord_y[0], coord_z[0], coord_x[3],
170  coord_y[3], coord_z[3], coord_x[2], coord_y[2],
171  coord_z[2], coord_px, coord_py, coord_pz,
172  &_dot_product[1], &_distance[1]);
173  calc_dot_product_tri(coord_x[0], coord_y[0], coord_z[0], coord_x[1],
174  coord_y[1], coord_z[1], coord_x[3], coord_y[3],
175  coord_z[3], coord_px, coord_py, coord_pz,
176  &_dot_product[2], &_distance[2]);
177  calc_dot_product_tri(coord_x[0], coord_y[0], coord_z[0], coord_x[2],
178  coord_y[2], coord_z[2], coord_x[1], coord_y[1],
179  coord_z[1], coord_px, coord_py, coord_pz,
180  &_dot_product[3], &_distance[3]);
181 
182  *dot_product = _dot_product[surf_id - 1];
183  *distance = _distance[surf_id - 1];
184 
185  n_positive = 0;
186  for (i = 0; i < 4; i++) {
187  if (_dot_product[i] > EPS_PRODUCT) {
188  if (_distance[i] > EPS_DISTANCE) {
189  return -1;
190  }
191  n_positive++;
192  }
193  }
194 
195  return n_positive;
196 }
197 
198 extern int HECMW_couple_judge_hex1(const struct hecmwST_local_mesh *local_mesh,
199  int elem, int surf_id, double coord_px,
200  double coord_py, double coord_pz,
201  double *dot_product, double *distance) {
202  double coord_x[8], coord_y[8], coord_z[8], _dot_product[8], _distance[8];
203  int node_index, node_id, n_positive, i;
204 
205  node_index = local_mesh->elem_node_index[elem - 1];
206  for (i = 0; i < 8; i++) {
207  node_id = local_mesh->elem_node_item[node_index + i];
208  coord_x[i] = local_mesh->node[3 * (node_id - 1)];
209  coord_y[i] = local_mesh->node[3 * (node_id - 1) + 1];
210  coord_z[i] = local_mesh->node[3 * (node_id - 1) + 2];
211  }
212 
213  calc_dot_product_quad(
214  coord_x[3], coord_y[3], coord_z[3], coord_x[0], coord_y[0], coord_z[0],
215  coord_x[4], coord_y[4], coord_z[4], coord_x[7], coord_y[7], coord_z[7],
216  coord_px, coord_py, coord_pz, &_dot_product[0], &_distance[0]);
217  calc_dot_product_quad(
218  coord_x[1], coord_y[1], coord_z[1], coord_x[2], coord_y[2], coord_z[2],
219  coord_x[6], coord_y[6], coord_z[6], coord_x[5], coord_y[5], coord_z[5],
220  coord_px, coord_py, coord_pz, &_dot_product[1], &_distance[1]);
221  calc_dot_product_quad(
222  coord_x[0], coord_y[0], coord_z[0], coord_x[1], coord_y[1], coord_z[1],
223  coord_x[5], coord_y[5], coord_z[5], coord_x[4], coord_y[4], coord_z[4],
224  coord_px, coord_py, coord_pz, &_dot_product[2], &_distance[2]);
225  calc_dot_product_quad(
226  coord_x[2], coord_y[2], coord_z[2], coord_x[3], coord_y[3], coord_z[3],
227  coord_x[7], coord_y[7], coord_z[7], coord_x[6], coord_y[6], coord_z[6],
228  coord_px, coord_py, coord_pz, &_dot_product[3], &_distance[3]);
229  calc_dot_product_quad(
230  coord_x[3], coord_y[3], coord_z[3], coord_x[2], coord_y[2], coord_z[2],
231  coord_x[1], coord_y[1], coord_z[1], coord_x[0], coord_y[0], coord_z[0],
232  coord_px, coord_py, coord_pz, &_dot_product[4], &_distance[4]);
233  calc_dot_product_quad(
234  coord_x[4], coord_y[4], coord_z[4], coord_x[5], coord_y[5], coord_z[5],
235  coord_x[6], coord_y[6], coord_z[6], coord_x[7], coord_y[7], coord_z[7],
236  coord_px, coord_py, coord_pz, &_dot_product[5], &_distance[5]);
237 
238  *dot_product = _dot_product[surf_id - 1];
239  *distance = _distance[surf_id - 1];
240 
241  n_positive = 0;
242  for (i = 0; i < 6; i++) {
243  if (_dot_product[i] > EPS_PRODUCT) {
244  if (_distance[i] > EPS_DISTANCE) {
245  return -1;
246  }
247  n_positive++;
248  }
249  }
250 
251  return n_positive;
252 }
hecmw_couple_judge.h
hecmwST_local_mesh::elem_node_item
int * elem_node_item
Definition: hecmw_struct.h:196
FRAC_1_3
#define FRAC_1_3
Definition: hecmw_couple_judge.c:21
hecmwST_local_mesh
Definition: hecmw_struct.h:139
HECMW_couple_judge_hex1
int HECMW_couple_judge_hex1(const struct hecmwST_local_mesh *local_mesh, int elem, int surf_id, double coord_px, double coord_py, double coord_pz, double *dot_product, double *distance)
Definition: hecmw_couple_judge.c:198
hecmw_struct.h
hecmwST_local_mesh::node
double * node
Definition: hecmw_struct.h:170
hecmwST_local_mesh::elem_node_index
int * elem_node_index
Definition: hecmw_struct.h:195
EPS_PRODUCT
#define EPS_PRODUCT
Definition: hecmw_couple_judge.c:25
hecmw_couple_struct.h
EPS_DISTANCE
#define EPS_DISTANCE
Definition: hecmw_couple_judge.c:26
FRAC_1_2
#define FRAC_1_2
Definition: hecmw_couple_judge.c:19
hecmw_couple_define.h
FRAC_1_4
#define FRAC_1_4
Definition: hecmw_couple_judge.c:23
HECMW_couple_judge_tet1
int HECMW_couple_judge_tet1(const struct hecmwST_local_mesh *local_mesh, int elem, int surf_id, double coord_px, double coord_py, double coord_pz, double *dot_product, double *distance)
Definition: hecmw_couple_judge.c:149