FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
30static 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
39static 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
49static 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
69static 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
87static 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
95static 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
121static 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
149extern 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
198extern 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}
#define FRAC_1_2
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)
#define EPS_PRODUCT
#define EPS_DISTANCE
#define FRAC_1_4
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)
#define FRAC_1_3