FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_couple_n2s_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 n2s_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 = node_id[i];
109 p->weight = 1.0 / (3.0 * area);
110 p->next = weight_list[i].next;
111 weight_list[id].next = p;
112 }
113
114 return 0;
115}
116
117static int n2s_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 = node_id[i];
146 p->weight = 1.0 / (4.0 * area);
147 p->next = weight_list[id].next;
148 weight_list[id].next = p;
149 }
150
151 return 0;
152}
153
154static int n2s_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->surf->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->surf->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
176 for (i = 0; i < boundary->surf->n; i++) {
177 elem = boundary->surf->item[2 * i];
178
179 if (mesh->elem_type[elem - 1] == HECMW_ETYPE_TET1) {
180 if (n2s_with_area_tet1(mesh, boundary, i, weight_list) != HECMW_SUCCESS)
181 goto error;
182 } else if (mesh->elem_type[elem - 1] == HECMW_ETYPE_HEX1) {
183 if (n2s_with_area_hex1(mesh, boundary, i, weight_list) != HECMW_SUCCESS)
184 goto error;
185 } else {
187 goto error;
188 }
189 }
190
191 /*
192 * make interpolating information
193 */
194 /* number of surfaces */
195 weight_info->n = boundary->surf->n;
196
197 /* interpolating type */
198 weight_info->type = HECMW_COUPLE_IP_NODE_TO_SURF;
199
200 /* index of list */
201 weight_info->index = (int *)HECMW_calloc(weight_info->n + 1, sizeof(int));
202 if (weight_info->index == NULL) {
203 HECMW_set_error(errno, "");
204 goto error;
205 }
206 for (n = 0, i = 0; i < boundary->surf->n; i++) {
207 for (p = weight_list[i].next; p; p = p->next) {
208 n++;
209 }
210 weight_info->index[i + 1] = n;
211 }
212
213 /* id & weight */
214 n_item = weight_info->index[weight_info->n];
215 weight_info->id = (int *)HECMW_malloc(sizeof(int) * n_item);
216 if (weight_info->id == NULL) {
217 HECMW_set_error(errno, "");
218 goto error;
219 }
220 weight_info->weight = (double *)HECMW_malloc(sizeof(double) * n_item);
221 if (weight_info->weight == NULL) {
222 HECMW_set_error(errno, "");
223 goto error;
224 }
225 for (n = 0, i = 0; i < boundary->surf->n; i++) {
226 for (p = weight_list[i].next; p; p = p->next) {
227 weight_info->id[n] = p->id;
228 weight_info->weight[n] = p->weight;
229 n++;
230 }
231 }
232
233 /*
234 * free linked list
235 */
236 for (i = 0; i < boundary->surf->n; i++) {
237 free_link_list(weight_list[i].next);
238 }
239 HECMW_free(weight_list);
240
241 return 0;
242
243error:
244 for (i = 0; i < boundary->surf->n; i++) {
245 free_link_list(weight_list[i].next);
246 }
247 HECMW_free(weight_list);
248 return -1;
249}
250
252 const struct hecmwST_local_mesh *mesh,
253 const struct hecmw_couple_boundary *boundary) {
254 struct hecmw_couple_weight_list *weight_info_list = NULL;
255 struct hecmw_couple_weight *weight_info = NULL;
256
257 if (mesh == NULL) {
259 "HECMW_couple_n2s_with_area(): 'mesh' is NULL");
260 return NULL;
261 }
262 if (boundary == NULL) {
264 "HECMW_couple_n2s_with_area(): 'boundary' is NULL");
265 return NULL;
266 }
267
268 if ((weight_info_list = HECMW_couple_alloc_weight_list()) == NULL)
269 return NULL;
270
271 if ((weight_info = HECMW_couple_alloc_weight()) == NULL) return NULL;
272 weight_info_list->info = weight_info;
273
274 if (n2s_with_area(mesh, boundary, weight_info)) goto error;
275
276 return weight_info_list;
277
278error:
279 HECMW_couple_free_weight(weight_info);
280 weight_info_list->info = NULL;
281 HECMW_couple_free_weight_list(weight_info_list);
282 return NULL;
283}
#define HECMW_ETYPE_HEX1
#define HECMW_ETYPE_TET1
#define HECMW_SUCCESS
Definition: hecmw_config.h:64
#define HECMW_COUPLE_IP_NODE_TO_SURF
#define HECMWCPL_E_NONSUPPORT_ETYPE
#define HECMWCPL_E_INVALID_ARG
#define FRAC_1_2
struct hecmw_couple_weight_list * HECMW_couple_n2s_with_area(const struct hecmwST_local_mesh *mesh, const struct hecmw_couple_boundary *boundary)
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