FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_couple_s2n_with_area.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"
18#include "hecmw_couple_weight.h"
20
21#define FRAC_1_2 (0.5)
22
23#define FRAC_1_3 (0.33333333333333333)
24
25#define FRAC_1_4 (0.25)
26
27/*================================================================================================*/
28
29struct link_list {
30 int id;
31 double weight;
32 struct link_list *next;
33};
34
36 double x;
37 double y;
38 double z;
39};
40
42 double x;
43 double y;
44 double z;
45};
46
47/*================================================================================================*/
48
49static void free_link_list(struct link_list *r) {
50 struct link_list *p, *q;
51
52 for (p = r; p; p = q) {
53 q = p->next;
54 HECMW_free(p);
55 }
56}
57
58static void cross_product(struct hecmw_couple_vertex *coord0,
59 struct hecmw_couple_vertex *coord1,
60 struct hecmw_couple_vertex *coord2,
61 struct hecmw_couple_vector *cross_prod) {
62 cross_prod->x = (coord1->y - coord0->y) * (coord2->z - coord0->z) -
63 (coord1->z - coord0->z) * (coord2->y - coord0->y);
64 cross_prod->y = (coord1->z - coord0->z) * (coord2->x - coord0->x) -
65 (coord1->x - coord0->x) * (coord2->z - coord0->z);
66 cross_prod->z = (coord1->x - coord0->x) * (coord2->y - coord0->y) -
67 (coord1->y - coord0->y) * (coord2->x - coord0->x);
68}
69
70static double tri_area(struct hecmw_couple_vector *cross_prod) {
71 return sqrt(cross_prod->x * cross_prod->x + cross_prod->y * cross_prod->y +
72 cross_prod->z * cross_prod->z) *
74}
75
76static double quad_area(struct hecmw_couple_vector *cross_prod1,
77 struct hecmw_couple_vector *cross_prod2) {
78 return tri_area(cross_prod1) + tri_area(cross_prod2);
79}
80
81static int s2n_with_area_tet1(const struct hecmwST_local_mesh *mesh,
82 const struct hecmw_couple_boundary *boundary,
83 int id, struct link_list *weight_list) {
84 struct hecmw_couple_vertex coord[3];
85 struct hecmw_couple_vector cross_prod;
86 struct link_list *p;
87 double area;
88 int node_id[3], node, n, i;
89
90 for (n = 0, i = boundary->elem_node_index[id];
91 i < boundary->elem_node_index[id + 1]; i++) {
92 node_id[n] = boundary->elem_node_item[i];
93 node = boundary->node->item[node_id[n]];
94 coord[n].x = mesh->node[3 * (node - 1)];
95 coord[n].y = mesh->node[3 * (node - 1) + 1];
96 coord[n].z = mesh->node[3 * (node - 1) + 2];
97 n++;
98 }
99 cross_product(&coord[0], &coord[1], &coord[2], &cross_prod);
100 area = tri_area(&cross_prod);
101
102 for (i = 0; i < 3; i++) {
103 p = (struct link_list *)HECMW_malloc(sizeof(struct link_list));
104 if (p == NULL) {
105 HECMW_set_error(errno, "");
106 return -1;
107 }
108 p->id = id;
109 p->weight = area * FRAC_1_3;
110 p->next = weight_list[node_id[i]].next;
111 weight_list[node_id[i]].next = p;
112 }
113
114 return 0;
115}
116
117static int s2n_with_area_hex1(const struct hecmwST_local_mesh *mesh,
118 const struct hecmw_couple_boundary *boundary,
119 int id, struct link_list *weight_list) {
120 struct hecmw_couple_vertex coord[4];
121 struct hecmw_couple_vector cross_prod1, cross_prod2;
122 struct link_list *p;
123 double area;
124 int node_id[4], node, n, i;
125
126 for (n = 0, i = boundary->elem_node_index[id];
127 i < boundary->elem_node_index[id + 1]; i++) {
128 node_id[n] = boundary->elem_node_item[i];
129 node = boundary->node->item[node_id[n]];
130 coord[n].x = mesh->node[3 * (node - 1)];
131 coord[n].y = mesh->node[3 * (node - 1) + 1];
132 coord[n].z = mesh->node[3 * (node - 1) + 2];
133 n++;
134 }
135 cross_product(&coord[0], &coord[1], &coord[3], &cross_prod1);
136 cross_product(&coord[2], &coord[3], &coord[1], &cross_prod2);
137 area = quad_area(&cross_prod1, &cross_prod2);
138
139 for (i = 0; i < 4; i++) {
140 p = (struct link_list *)HECMW_malloc(sizeof(struct link_list));
141 if (p == NULL) {
142 HECMW_set_error(errno, "");
143 return -1;
144 }
145 p->id = id;
146 p->weight = area * FRAC_1_4;
147 p->next = weight_list[node_id[i]].next;
148 weight_list[node_id[i]].next = p;
149 }
150
151 return 0;
152}
153
154static int s2n_with_area(const struct hecmwST_local_mesh *mesh,
155 const struct hecmw_couple_boundary *boundary,
156 struct hecmw_couple_weight *weight_info) {
157 struct link_list *weight_list = NULL, *p;
158 int elem, n_item, size, n, i;
159
160 size = sizeof(struct link_list) * boundary->node->n;
161 weight_list = (struct link_list *)HECMW_malloc(size);
162 if (weight_list == NULL) {
163 HECMW_set_error(errno, "");
164 goto error;
165 }
166 for (i = 0; i < boundary->node->n; i++) {
167 weight_list[i].id = -1;
168 weight_list[i].weight = 0.0;
169 weight_list[i].next = NULL;
170 }
171
172 /*
173 * calculate weight
174 */
175 for (i = 0; i < boundary->surf->n; i++) {
176 elem = boundary->surf->item[2 * i];
177
178 if (mesh->elem_type[elem - 1] == HECMW_ETYPE_TET1) {
179 if (s2n_with_area_tet1(mesh, boundary, i, weight_list)) goto error;
180 } else if (mesh->elem_type[elem - 1] == HECMW_ETYPE_HEX1) {
181 if (s2n_with_area_hex1(mesh, boundary, i, weight_list)) goto error;
182 } else {
184 goto error;
185 }
186 }
187
188 /*
189 * make interpolating information
190 */
191 /* number of nodes */
192 weight_info->n = boundary->node->n;
193
194 /* interpolating type */
195 weight_info->type = HECMW_COUPLE_IP_SURF_TO_NODE;
196
197 /* index of list */
198 weight_info->index = (int *)HECMW_calloc(weight_info->n + 1, sizeof(int));
199 if (weight_info->index == NULL) {
200 HECMW_set_error(errno, "");
201 goto error;
202 }
203 for (n = 0, i = 0; i < boundary->node->n; i++) {
204 for (p = weight_list[i].next; p; p = p->next) {
205 n++;
206 }
207 weight_info->index[i + 1] = n;
208 }
209
210 /* id which interpolates nodes and its weight */
211 n_item = weight_info->index[weight_info->n];
212 weight_info->id = (int *)HECMW_malloc(sizeof(int) * n_item);
213 if (weight_info->id == NULL) {
214 HECMW_set_error(errno, "");
215 goto error;
216 }
217 weight_info->weight = (double *)HECMW_malloc(sizeof(double) * n_item);
218 if (weight_info->weight == NULL) {
219 HECMW_set_error(errno, "");
220 goto error;
221 }
222 for (n = 0, i = 0; i < boundary->node->n; i++) {
223 for (p = weight_list[i].next; p; p = p->next) {
224 weight_info->id[n] = p->id;
225 weight_info->weight[n] = p->weight;
226 n++;
227 }
228 }
229
230 /*
231 * free linked list
232 */
233 for (i = 0; i < boundary->node->n; i++) {
234 free_link_list(weight_list[i].next);
235 }
236 HECMW_free(weight_list);
237
238 return 0;
239
240error:
241 for (i = 0; i < boundary->node->n; i++) {
242 free_link_list(weight_list[i].next);
243 }
244 HECMW_free(weight_list);
245 return -1;
246}
247
249 const struct hecmwST_local_mesh *mesh,
250 const struct hecmw_couple_boundary *boundary) {
251 struct hecmw_couple_weight_list *weight_info_list = NULL;
252 struct hecmw_couple_weight *weight_info = NULL;
253 struct link_list *weight_list = NULL, *p;
254 int elem, n_item, size, n, i;
255
256 if (mesh == NULL) {
258 "HECMW_couple_s2n_by_area(): 'mesh' is NULL");
259 return NULL;
260 }
261 if (boundary == NULL) {
263 "HECMW_couple_s2n_by_area(): 'boundary' is NULL");
264 return NULL;
265 }
266
267 if ((weight_info_list = HECMW_couple_alloc_weight_list()) == NULL)
268 return NULL;
269
270 if ((weight_info = HECMW_couple_alloc_weight()) == NULL) goto error;
271 weight_info_list->info = weight_info;
272
273 if (s2n_with_area(mesh, boundary, weight_info)) goto error;
274
275 return weight_info_list;
276
277error:
278 HECMW_couple_free_weight(weight_info);
279 weight_info_list->info = NULL;
280 HECMW_couple_free_weight_list(weight_info_list);
281 return NULL;
282}
#define HECMW_ETYPE_HEX1
#define HECMW_ETYPE_TET1
#define HECMW_COUPLE_IP_SURF_TO_NODE
#define HECMWCPL_E_NONSUPPORT_ETYPE
#define HECMWCPL_E_INVALID_ARG
struct hecmw_couple_weight_list * HECMW_couple_s2n_with_area(const struct hecmwST_local_mesh *mesh, const struct hecmw_couple_boundary *boundary)
#define FRAC_1_2
#define FRAC_1_4
#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)
struct hecmwST_local_mesh * mesh
Definition: hecmw_repart.h:71
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
subroutine cross_product(v1, v2, vn)
Definition: utilities.f90:330
struct hecmw_couple_boundary_item * node
struct hecmw_couple_boundary_item * surf
struct hecmw_couple_weight * info