FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_couple_s2n_dist_surf.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 <errno.h>
10#include <assert.h>
11
12#include "hecmw_common_define.h"
13#include "hecmw_struct.h"
14#include "hecmw_malloc.h"
15
16#include "hecmw_couple_define.h"
17#include "hecmw_couple_struct.h"
18#include "hecmw_couple_weight.h"
23
24#define FRAC_1_2 (0.5)
25
26#define FRAC_1_3 (0.33333333333333333)
27
28#define FRAC_1_4 (0.25)
29
30#define EPS_ZERO (1.0E-24)
31
32struct link_list {
33 int id;
34 double weight;
35 struct link_list *next;
36};
37
39 double x;
40 double y;
41 double z;
42};
43
45 double x;
46 double y;
47 double z;
48};
49
50/*================================================================================================*/
51
52static void free_link_list(struct link_list *r) {
53 struct link_list *p, *q;
54
55 for (p = r; p; p = q) {
56 q = p->next;
57 HECMW_free(p);
58 }
59}
60
61static int intercomm_d2s_coord(
62 const struct hecmw_couple_mapped_point *mapped_point,
63 const struct hecmw_couple_inter_iftable *inter_tbl,
64 const struct hecmw_couple_comm *comm_src,
65 const struct hecmw_couple_comm *comm_dst,
66 const struct hecmw_couple_comm *intercomm, double **coord) {
67 int *sendbuf_index = NULL, *recvbuf_index = NULL;
68 double *sendbuf = NULL;
69 int size, rtc, i;
70
71 if (comm_dst->is_member) {
72 sendbuf_index =
73 (int *)HECMW_calloc(inter_tbl->n_neighbor_pe_import + 1, sizeof(int));
74 if (sendbuf_index == NULL) {
75 HECMW_set_error(errno, "");
76 goto error;
77 }
78 for (i = 0; i <= inter_tbl->n_neighbor_pe_import; i++) {
79 sendbuf_index[i] = 3 * inter_tbl->import_index[i];
80 }
81
82 size = sizeof(double) *
83 inter_tbl->import_index[inter_tbl->n_neighbor_pe_import] * 3 +
84 1;
85 sendbuf = (double *)HECMW_malloc(size);
86 if (sendbuf == NULL) {
87 HECMW_set_error(errno, "");
88 goto error;
89 }
90 for (i = 0; i < inter_tbl->import_index[inter_tbl->n_neighbor_pe_import];
91 i++) {
92 sendbuf[3 * i] = mapped_point->coord[3 * (inter_tbl->import_item[i])];
93 sendbuf[3 * i + 1] =
94 mapped_point->coord[3 * (inter_tbl->import_item[i]) + 1];
95 sendbuf[3 * i + 2] =
96 mapped_point->coord[3 * (inter_tbl->import_item[i]) + 2];
97 }
98 }
99
100 if (comm_src->is_member) {
101 recvbuf_index =
102 (int *)HECMW_calloc(inter_tbl->n_neighbor_pe_export + 1, sizeof(int));
103 if (recvbuf_index == NULL) {
104 HECMW_set_error(errno, "");
105 goto error;
106 }
107 for (i = 0; i <= inter_tbl->n_neighbor_pe_export; i++) {
108 recvbuf_index[i] = 3 * inter_tbl->export_index[i];
109 }
110
111 size = sizeof(double) *
112 inter_tbl->export_index[inter_tbl->n_neighbor_pe_export] * 3 +
113 1;
114 *coord = (double *)HECMW_malloc(size);
115 if (*coord == NULL) {
116 HECMW_set_error(errno, "");
117 goto error;
118 }
119 }
120
122 inter_tbl->n_neighbor_pe_import, inter_tbl->neighbor_pe_import,
123 sendbuf_index, sendbuf, inter_tbl->n_neighbor_pe_export,
124 inter_tbl->neighbor_pe_export, recvbuf_index, *coord, HECMW_DOUBLE,
125 intercomm->comm);
126 if (rtc != 0) goto error;
127
128 HECMW_free(sendbuf_index);
129 HECMW_free(sendbuf);
130 HECMW_free(recvbuf_index);
131 return 0;
132
133error:
134 HECMW_free(sendbuf_index);
135 HECMW_free(sendbuf);
136 HECMW_free(recvbuf_index);
137 return -1;
138}
139
140static int s2n_dist_surf_tet1(const struct hecmwST_local_mesh *mesh_src,
141 const struct hecmw_couple_boundary *boundary_src,
142 int id,
143 const struct hecmw_couple_vertex *coord_dst,
144 struct link_list *weight_list_surf,
145 struct link_list *weight_list_node) {
146 struct link_list *p;
147 struct hecmw_couple_vertex coord[3], gravity;
148 double d, r_d_surf, r_d_node[3], r_d_sum = 0.0;
149 int node_id[3], node, n, i;
150
151 for (n = 0, i = boundary_src->elem_node_index[id];
152 i < boundary_src->elem_node_index[id + 1]; i++) {
153 node_id[n] = boundary_src->elem_node_item[i];
154 node = boundary_src->node->item[node_id[n]];
155 coord[n].x = mesh_src->node[3 * (node - 1)];
156 coord[n].y = mesh_src->node[3 * (node - 1) + 1];
157 coord[n].z = mesh_src->node[3 * (node - 1) + 2];
158 n++;
159 }
160 gravity.x = (coord[0].x + coord[1].x + coord[2].x) * FRAC_1_3;
161 gravity.y = (coord[0].y + coord[1].y + coord[2].y) * FRAC_1_3;
162 gravity.z = (coord[0].z + coord[1].z + coord[2].z) * FRAC_1_3;
163
164 d = sqrt((gravity.x - coord_dst->x) * (gravity.x - coord_dst->x) +
165 (gravity.y - coord_dst->y) * (gravity.y - coord_dst->y) +
166 (gravity.z - coord_dst->z) * (gravity.z - coord_dst->z));
167 r_d_surf = 1.0 / (d + EPS_ZERO);
168 r_d_sum += r_d_surf;
169
170 for (i = 0; i < 3; i++) {
171 d = sqrt((coord[i].x - coord_dst->x) * (coord[i].x - coord_dst->x) +
172 (coord[i].y - coord_dst->y) * (coord[i].y - coord_dst->y) +
173 (coord[i].z - coord_dst->z) * (coord[i].z - coord_dst->z));
174 r_d_node[i] = 1.0 / (d + EPS_ZERO);
175 r_d_sum += r_d_node[i];
176 }
177
178 p = (struct link_list *)HECMW_malloc(sizeof(struct link_list));
179 if (p == NULL) {
180 HECMW_set_error(errno, "");
181 return -1;
182 }
183 p->id = id;
184 p->weight = r_d_surf / r_d_sum;
185 p->next = weight_list_surf->next;
186 weight_list_surf->next = p;
187
188 for (i = 0; i < 3; i++) {
189 p = (struct link_list *)HECMW_malloc(sizeof(struct link_list));
190 if (p == NULL) {
191 HECMW_set_error(errno, "");
192 return -1;
193 }
194 p->id = node_id[i];
195 p->weight = r_d_node[i] / r_d_sum;
196 p->next = weight_list_node->next;
197 weight_list_node->next = p;
198 }
199
200 return 0;
201}
202
203static int s2n_dist_surf_hex1(const struct hecmwST_local_mesh *mesh_src,
204 const struct hecmw_couple_boundary *boundary_src,
205 int id,
206 const struct hecmw_couple_vertex *coord_dst,
207 struct link_list *weight_list_surf,
208 struct link_list *weight_list_node) {
209 struct link_list *p;
210 struct hecmw_couple_vertex coord[4], gravity;
211 double d, r_d_surf, r_d_node[4], r_d_sum = 0.0;
212 int node_id[4], node, n, i;
213
214 for (n = 0, i = boundary_src->elem_node_index[id];
215 i < boundary_src->elem_node_index[id + 1]; i++) {
216 node_id[n] = boundary_src->elem_node_item[i];
217 node = boundary_src->node->item[node_id[n]];
218 coord[n].x = mesh_src->node[3 * (node - 1)];
219 coord[n].y = mesh_src->node[3 * (node - 1) + 1];
220 coord[n].z = mesh_src->node[3 * (node - 1) + 2];
221 n++;
222 }
223 gravity.x = (coord[0].x + coord[1].x + coord[2].x + coord[3].x) * FRAC_1_4;
224 gravity.y = (coord[0].y + coord[1].y + coord[2].y + coord[3].y) * FRAC_1_4;
225 gravity.z = (coord[0].z + coord[1].z + coord[2].z + coord[3].z) * FRAC_1_4;
226
227 d = sqrt((gravity.x - coord_dst->x) * (gravity.x - coord_dst->x) +
228 (gravity.y - coord_dst->y) * (gravity.y - coord_dst->y) +
229 (gravity.z - coord_dst->z) * (gravity.z - coord_dst->z));
230 r_d_surf = 1.0 / (d + EPS_ZERO);
231 r_d_sum += r_d_surf;
232
233 for (i = 0; i < 4; i++) {
234 d = sqrt((coord[i].x - coord_dst->x) * (coord[i].x - coord_dst->x) +
235 (coord[i].y - coord_dst->y) * (coord[i].y - coord_dst->y) +
236 (coord[i].z - coord_dst->z) * (coord[i].z - coord_dst->z));
237 r_d_node[i] = 1.0 / (d + EPS_ZERO);
238 r_d_sum += r_d_node[i];
239 }
240
241 p = (struct link_list *)HECMW_malloc(sizeof(struct link_list));
242 if (p == NULL) {
243 HECMW_set_error(errno, "");
244 return -1;
245 }
246 p->id = id;
247 p->weight = r_d_surf / r_d_sum;
248 p->next = weight_list_surf->next;
249 weight_list_surf->next = p;
250
251 for (i = 0; i < 4; i++) {
252 p = (struct link_list *)HECMW_malloc(sizeof(struct link_list));
253 if (p == NULL) {
254 HECMW_set_error(errno, "");
255 return -1;
256 }
257 p->id = node_id[i];
258 p->weight = r_d_node[i] / r_d_sum;
259 p->next = weight_list_node->next;
260 weight_list_node->next = p;
261 }
262
263 return 0;
264}
265
266static int s2n_dist_surf(const struct hecmwST_local_mesh *mesh_src,
267 const struct hecmw_couple_boundary *boundary_src,
268 const struct hecmw_couple_inter_iftable *inter_tbl,
269 double *coord,
270 struct hecmw_couple_weight *weight_info_surf,
271 struct hecmw_couple_weight *weight_info_node) {
272 struct link_list *weight_list_surf = NULL, *weight_list_node = NULL, *p;
273 struct hecmw_couple_vertex coord_dst;
274 int elem, n_item, id, rtc, n, i;
275
276 n_item = inter_tbl->export_index[inter_tbl->n_neighbor_pe_export] + 1;
277 weight_list_surf =
278 (struct link_list *)HECMW_malloc(sizeof(struct link_list) * n_item);
279 if (weight_list_surf == NULL) {
280 HECMW_set_error(errno, "");
281 goto error;
282 }
283 weight_list_node =
284 (struct link_list *)HECMW_malloc(sizeof(struct link_list) * n_item);
285 if (weight_list_node == NULL) {
286 HECMW_set_error(errno, "");
287 goto error;
288 }
289 for (i = 0; i < n_item; i++) {
290 weight_list_surf[i].id = -1;
291 weight_list_surf[i].weight = 0.0;
292 weight_list_surf[i].next = NULL;
293 weight_list_node[i].id = -1;
294 weight_list_node[i].weight = 0.0;
295 weight_list_node[i].next = NULL;
296 }
297
298 for (i = 0; i < inter_tbl->export_index[inter_tbl->n_neighbor_pe_export];
299 i++) {
300 coord_dst.x = coord[3 * i];
301 coord_dst.y = coord[3 * i + 1];
302 coord_dst.z = coord[3 * i + 2];
303
304 id = inter_tbl->export_item[i];
305 elem = boundary_src->surf->item[2 * id];
306 if (mesh_src->elem_type[elem - 1] == HECMW_ETYPE_TET1) {
307 rtc = s2n_dist_surf_tet1(mesh_src, boundary_src, id, &coord_dst,
308 &weight_list_surf[i], &weight_list_node[i]);
309 if (rtc != HECMW_SUCCESS) goto error;
310 } else if (mesh_src->elem_type[elem - 1] == HECMW_ETYPE_HEX1) {
311 rtc = s2n_dist_surf_hex1(mesh_src, boundary_src, id, &coord_dst,
312 &weight_list_surf[i], &weight_list_node[i]);
313 if (rtc != HECMW_SUCCESS) goto error;
314 } else {
316 goto error;
317 }
318 }
319
320 weight_info_node->n =
321 inter_tbl->export_index[inter_tbl->n_neighbor_pe_export];
322 weight_info_node->type = HECMW_COUPLE_IP_NODE_TO_NODE;
323
324 weight_info_node->index =
325 (int *)HECMW_calloc(weight_info_node->n + 1, sizeof(int));
326 if (weight_info_node->index == NULL) {
327 HECMW_set_error(errno, "");
328 goto error;
329 }
330 for (n = 0, i = 0;
331 i < inter_tbl->export_index[inter_tbl->n_neighbor_pe_export]; i++) {
332 for (p = weight_list_node[i].next; p; p = p->next) {
333 n++;
334 }
335 weight_info_node->index[i + 1] = n;
336 }
337
338 n_item = weight_info_node->index[weight_info_node->n];
339 weight_info_node->id = (int *)HECMW_malloc(sizeof(int) * n_item);
340 if (weight_info_node->id == NULL) {
341 HECMW_set_error(errno, "");
342 goto error;
343 }
344 weight_info_node->weight = (double *)HECMW_malloc(sizeof(double) * n_item);
345 if (weight_info_node->weight == NULL) {
346 HECMW_set_error(errno, "");
347 goto error;
348 }
349 for (n = 0, i = 0;
350 i < inter_tbl->export_index[inter_tbl->n_neighbor_pe_export]; i++) {
351 for (p = weight_list_node[i].next; p; p = p->next) {
352 weight_info_node->id[n] = p->id;
353 weight_info_node->weight[n] = p->weight;
354 n++;
355 }
356 }
357
358 weight_info_surf->n =
359 inter_tbl->export_index[inter_tbl->n_neighbor_pe_export];
360 weight_info_surf->type = HECMW_COUPLE_IP_SURF_TO_NODE;
361
362 weight_info_surf->index =
363 (int *)HECMW_calloc(weight_info_surf->n + 1, sizeof(int));
364 if (weight_info_surf->index == NULL) {
365 HECMW_set_error(errno, "");
366 goto error;
367 }
368 for (n = 0, i = 0;
369 i < inter_tbl->export_index[inter_tbl->n_neighbor_pe_export]; i++) {
370 for (p = weight_list_surf[i].next; p; p = p->next) {
371 n++;
372 }
373 weight_info_surf->index[i + 1] = n;
374 }
375
376 n_item = weight_info_surf->index[weight_info_node->n];
377 weight_info_surf->id = (int *)HECMW_malloc(sizeof(int) * n_item);
378 if (weight_info_surf->id == NULL) {
379 HECMW_set_error(errno, "");
380 goto error;
381 }
382 weight_info_surf->weight = (double *)HECMW_malloc(sizeof(double) * n_item);
383 if (weight_info_surf->weight == NULL) {
384 HECMW_set_error(errno, "");
385 goto error;
386 }
387 for (n = 0, i = 0;
388 i < inter_tbl->export_index[inter_tbl->n_neighbor_pe_export]; i++) {
389 for (p = weight_list_surf[i].next; p; p = p->next) {
390 weight_info_surf->id[n] = p->id;
391 weight_info_surf->weight[n] = p->weight;
392 n++;
393 }
394 }
395
396 for (i = 0; i < inter_tbl->export_index[inter_tbl->n_neighbor_pe_export];
397 i++) {
398 free_link_list(weight_list_surf[i].next);
399 free_link_list(weight_list_node[i].next);
400 }
401 HECMW_free(weight_list_surf);
402 HECMW_free(weight_list_node);
403
404 return 0;
405
406error:
407 for (i = 0; i < inter_tbl->export_index[inter_tbl->n_neighbor_pe_export];
408 i++) {
409 free_link_list(weight_list_surf[i].next);
410 free_link_list(weight_list_node[i].next);
411 }
412 HECMW_free(weight_list_surf);
413 HECMW_free(weight_list_node);
414
415 return -1;
416}
417
419 const struct hecmwST_local_mesh *mesh_src,
420 const struct hecmwST_local_mesh *mesh_dst,
421 const struct hecmw_couple_comm *comm_src,
422 const struct hecmw_couple_comm *comm_dst,
423 const struct hecmw_couple_comm *intercomm,
424 const struct hecmw_couple_boundary *boundary_src,
425 const struct hecmw_couple_boundary *boundary_dst,
426 const struct hecmw_couple_mapped_point *mapped_point,
427 const struct hecmw_couple_inter_iftable *inter_tbl) {
428 struct hecmw_couple_weight_list *weight_info_list = NULL;
429 struct hecmw_couple_weight *weight_info_node = NULL;
430 struct hecmw_couple_weight *weight_info_surf = NULL;
431 double *coord = NULL;
432 int rtc, i;
433
434 if ((weight_info_list = HECMW_couple_alloc_weight_list()) == NULL)
435 return NULL;
436
437 rtc = intercomm_d2s_coord(mapped_point, inter_tbl, comm_src, comm_dst,
438 intercomm, &coord);
439 if (rtc) goto error;
440
441 if (comm_src->is_member) {
442 if ((weight_info_list->next = HECMW_couple_alloc_weight_list()) == NULL)
443 goto error;
444
445 if ((weight_info_node = HECMW_couple_alloc_weight()) == NULL) goto error;
446 if ((weight_info_surf = HECMW_couple_alloc_weight()) == NULL) goto error;
447
448 weight_info_list->info = weight_info_node;
449 weight_info_list->next->info = weight_info_surf;
450
451 rtc = s2n_dist_surf(mesh_src, boundary_src, inter_tbl, coord,
452 weight_info_surf, weight_info_node);
453 if (rtc) goto error;
454 }
455
456 return weight_info_list;
457
458error:
459 HECMW_couple_free_weight(weight_info_node);
460 HECMW_couple_free_weight(weight_info_surf);
461 HECMW_couple_free_weight_list(weight_info_list);
462 return NULL;
463}
#define HECMW_ETYPE_HEX1
#define HECMW_ETYPE_TET1
#define HECMW_DOUBLE
Definition: hecmw_config.h:50
#define HECMW_SUCCESS
Definition: hecmw_config.h:64
int HECMW_couple_inter_send_recv(int n_neighbor_pe_send, int *neighbor_pe_send, int *sendbuf_index, void *sendbuf, int n_neighbor_pe_recv, int *neighbor_pe_recv, int *recvbuf_index, void *recvbuf, HECMW_Datatype datatype, HECMW_Comm comm)
#define HECMW_COUPLE_IP_NODE_TO_NODE
#define HECMW_COUPLE_IP_SURF_TO_NODE
#define HECMWCPL_E_NONSUPPORT_ETYPE
#define EPS_ZERO
#define FRAC_1_4
struct hecmw_couple_weight_list * HECMW_couple_s2n_dist_surf(const struct hecmwST_local_mesh *mesh_src, const struct hecmwST_local_mesh *mesh_dst, const struct hecmw_couple_comm *comm_src, const struct hecmw_couple_comm *comm_dst, const struct hecmw_couple_comm *intercomm, const struct hecmw_couple_boundary *boundary_src, const struct hecmw_couple_boundary *boundary_dst, const struct hecmw_couple_mapped_point *mapped_point, const struct hecmw_couple_inter_iftable *inter_tbl)
#define FRAC_1_3
struct hecmw_couple_weight * HECMW_couple_alloc_weight(void)
struct hecmw_couple_weight_list * HECMW_couple_alloc_weight_list(void)
void HECMW_couple_free_weight(struct hecmw_couple_weight *p)
void HECMW_couple_free_weight_list(struct hecmw_couple_weight_list *r)
int HECMW_set_error(int errorno, const char *fmt,...)
Definition: hecmw_error.c:37
#define NULL
#define HECMW_calloc(nmemb, size)
Definition: hecmw_malloc.h:21
#define HECMW_free(ptr)
Definition: hecmw_malloc.h:24
#define HECMW_malloc(size)
Definition: hecmw_malloc.h:20
struct hecmw_couple_boundary_item * node
struct hecmw_couple_boundary_item * surf
struct hecmw_couple_weight_list * next
struct hecmw_couple_weight * info