1 /*
2 * authors: Claus-Justus Heine
3 * Abteilung fuer Angewandte Mathematik
4 * Albert-Ludwigs-Universit"at Freiburg
5 * Hermann-Herder-Str. 10
6 * 79104 Freiburg
7 * Germany
8 * claus@mathematik.uni-freiburg.de
9 *
10 * Copyright (c) by C.-J. Heine (2005-2007)
11 *
12 * Generic support for assembling linear systems for time dependent problems.
13 */
14
15 #include "alberta_intern.h"
16 #include "alberta.h"
17
18 /** Theta-splitting scheme, given the element mass matrix and the part
19 * of the stiffness matrix which is subject to the \theta-splitting.
20 *
21 * Parameters:
22 *
23 * @param[in] el_info Element descriptor, filled as specified by
24 * ::EL_MATRIX_INFO::fill_flag.
25 * @param[in,out] dof_matrix The (global) system matrix.
26 * @param[in,out] f_h The righ-hand-side vector.
27 * @param[in] u_h The interpolation of the old solution onto the new
28 * mesh.
29 * @param[in] tau Time-step size.
30 * @param[in] mass_info The fill-info for the element mass matrix.
31 * @param[in] theta The splitting coefficient for the elliptic part.
32 * @param[in] stiff_info The fill-info for the element stiffness matrix,
33 * the elliptic part of the PDE.
34 *
35 */
36
37 typedef struct system_fill_info_instat
38 {
39 EL_SYS_INFO_INSTAT elsii;
40 const DOF_REAL_VEC *u_h;
41 const EL_MATRIX_INFO *stiff_info;
42 const EL_MATRIX_INFO *mass_info;
43 int n_row;
44 int n_col;
45 int n_row_max;
46 int n_col_max;
47 } SYSTEM_FILL_INFO_INSTAT;
48
49 typedef struct system_fill_info_dow_instat
50 {
51 EL_SYS_INFO_DOW_INSTAT elsii;
52 const DOF_REAL_VEC_D *u_h_d;
53 const EL_MATRIX_INFO *stiff_info;
54 const EL_MATRIX_INFO *mass_info;
55 int n_row;
56 int n_col;
57 int n_row_max;
58 int n_col_max;
59 } SYSTEM_FILL_INFO_DOW_INSTAT;
60
element_system_instat(const EL_INFO * el_info,REAL tau,REAL theta,EL_SYS_INFO_INSTAT * thisptr)61 static INIT_EL_TAG element_system_instat(const EL_INFO *el_info,
62 REAL tau,
63 REAL theta,
64 EL_SYS_INFO_INSTAT *thisptr)
65 {
66 SYSTEM_FILL_INFO_INSTAT *info = (SYSTEM_FILL_INFO_INSTAT *)thisptr;
67 REAL theta1 = 1.0 - theta;
68 REAL tau1 = 1.0/tau;
69 int n_row, n_col;
70
71 /* First compute the element matrices, el_matrix_fct() will call all
72 * necessary INIT_ELEMENT() methods for us.
73 */
74 thisptr->el_stiff =
75 info->stiff_info->el_matrix_fct(el_info, info->stiff_info->fill_info);
76 if (thisptr->el_stiff == NULL) {
77 return INIT_EL_TAG_NULL;
78 }
79 thisptr->el_mass =
80 info->mass_info->el_matrix_fct(el_info, info->mass_info->fill_info);
81 if (thisptr->el_mass == NULL) {
82 return INIT_EL_TAG_NULL;
83 }
84
85 thisptr->el_matrix->n_row =
86 n_row = thisptr->row_fe_space->bas_fcts->n_bas_fcts;
87 thisptr->el_matrix->n_col =
88 n_col = thisptr->col_fe_space->bas_fcts->n_bas_fcts;
89
90 info->n_row = n_row;
91 info->n_col = n_col;
92
93 fill_el_real_vec((EL_REAL_VEC *)thisptr->u_h_loc, el_info->el, info->u_h);
94
95 /* f += -(1-theta)*a(u_old,psi_i) + 1/tau*m(u_old,psi_i) */
96 if (theta1) {
97 el_bi_mat_vec(-theta1, thisptr->el_stiff, tau1, thisptr->el_mass,
98 thisptr->u_h_loc, 1.0, thisptr->el_load);
99 } else {
100 el_mat_vec(tau1, thisptr->el_mass, thisptr->u_h_loc, thisptr->el_load);
101 }
102
103 /* initialize with 1/tau*m(psi_i,psi_j) */
104 el_mat_axey(tau1, thisptr->el_mass, thisptr->el_matrix);
105
106 /* add theta*a(psi_i,psi_j) */
107 if (theta) {
108 el_mat_axpy(theta, thisptr->el_stiff, thisptr->el_matrix);
109 }
110 return INIT_EL_TAG_DFLT;
111 }
112
element_system_dow_instat(const EL_INFO * el_info,REAL tau,REAL theta,EL_SYS_INFO_DOW_INSTAT * thisptr)113 static INIT_EL_TAG element_system_dow_instat(const EL_INFO *el_info,
114 REAL tau,
115 REAL theta,
116 EL_SYS_INFO_DOW_INSTAT *thisptr)
117 {
118 SYSTEM_FILL_INFO_DOW_INSTAT *info = (SYSTEM_FILL_INFO_DOW_INSTAT *)thisptr;
119 REAL theta1 = 1.0 - theta;
120 REAL tau1 = 1.0/tau;
121 int n_row, n_col;
122
123 /* First compute the element matrices, el_matrix_fct() will call all
124 * necessary INIT_ELEMENT() methods for us.
125 */
126 thisptr->el_stiff = info->stiff_info->el_matrix_fct(el_info,
127 info->stiff_info->fill_info);
128 if (thisptr->el_stiff == NULL) {
129 return INIT_EL_TAG_NULL;
130 }
131
132 thisptr->el_mass =
133 info->mass_info->el_matrix_fct(el_info, info->mass_info->fill_info);
134 if (thisptr->el_mass == NULL) {
135 return INIT_EL_TAG_NULL;
136 }
137
138 thisptr->el_matrix->n_row =
139 n_row = thisptr->row_fe_space->bas_fcts->n_bas_fcts;
140 thisptr->el_matrix->n_col =
141 n_col = thisptr->col_fe_space->bas_fcts->n_bas_fcts;
142
143 info->n_row = n_row;
144 info->n_col = n_col;
145
146 fill_el_real_vec_d(
147 (EL_REAL_VEC_D *)thisptr->u_h_loc, el_info->el, info->u_h_d);
148
149 /* f += -(1-theta)*a(u_old,psi_i) + 1/tau*m(u_old,psi_i) */
150 if (theta1) {
151 el_bi_mat_vec_dow(-theta1, thisptr->el_stiff, tau1, thisptr->el_mass,
152 thisptr->u_h_loc, 1.0, thisptr->el_load);
153 } else {
154 el_mat_vec_dow(tau1, thisptr->el_mass, thisptr->u_h_loc, thisptr->el_load);
155 }
156
157 /* initialize with 1/tau*m(psi_i,psi_j) */
158 el_mat_axey(tau1, thisptr->el_mass, thisptr->el_matrix);
159
160 /* add theta*a(psi_i,psi_j) */
161 if (theta) {
162 el_mat_axpy(theta, thisptr->el_stiff, thisptr->el_matrix);
163 }
164 return INIT_EL_TAG_DFLT;
165 }
166
fill_sys_info_instat(const OPERATOR_INFO * stiff_info,const OPERATOR_INFO * mass_info,const DOF_REAL_VEC * u_h)167 EL_SYS_INFO_INSTAT *fill_sys_info_instat(const OPERATOR_INFO *stiff_info,
168 const OPERATOR_INFO *mass_info,
169 const DOF_REAL_VEC *u_h)
170 {
171 SYSTEM_FILL_INFO_INSTAT *info;
172
173 info = MEM_CALLOC(1, SYSTEM_FILL_INFO_INSTAT);
174
175 info->stiff_info = fill_matrix_info(stiff_info, NULL);
176 info->mass_info = fill_matrix_info(mass_info, NULL);
177
178 info->elsii.row_fe_space = info->mass_info->row_fe_space;
179 info->elsii.col_fe_space = info->mass_info->col_fe_space;
180 if (info->elsii.col_fe_space == NULL) {
181 info->elsii.col_fe_space = info->elsii.row_fe_space;
182 }
183
184 info->elsii.el_update_fct = element_system_instat;
185
186 info->n_row = info->elsii.row_fe_space->bas_fcts->n_bas_fcts;
187 info->n_row_max = info->elsii.row_fe_space->bas_fcts->n_bas_fcts_max;
188 info->n_col = info->elsii.col_fe_space->bas_fcts->n_bas_fcts;
189 info->n_col_max = info->elsii.col_fe_space->bas_fcts->n_bas_fcts_max;
190
191 info->elsii.el_matrix = get_el_matrix(info->elsii.row_fe_space,
192 info->elsii.col_fe_space,
193 MATENT_REAL);
194 info->elsii.el_load = get_el_real_vec(info->elsii.row_fe_space->bas_fcts);
195 info->elsii.u_h_loc = get_el_real_vec(info->elsii.col_fe_space->bas_fcts);
196
197 info->elsii.fill_flag =
198 info->stiff_info->fill_flag | info->mass_info->fill_flag;
199 BNDRY_FLAGS_CPY(info->elsii.dirichlet_bndry,
200 info->mass_info->dirichlet_bndry);
201 BNDRY_FLAGS_OR(info->elsii.dirichlet_bndry,
202 info->stiff_info->dirichlet_bndry);
203 if (!BNDRY_FLAGS_IS_INTERIOR(info->elsii.dirichlet_bndry)) {
204 info->elsii.fill_flag |= FILL_BOUND;
205 if (info->elsii.row_fe_space->mesh->is_periodic &&
206 !(info->elsii.row_fe_space->admin->flags & ADM_PERIODIC))
207 info->elsii.fill_flag |= FILL_NON_PERIODIC;
208 }
209
210 info->u_h = u_h;
211
212 return &info->elsii;
213 }
214
215 EL_SYS_INFO_DOW_INSTAT *
fill_sys_info_instat_dow(const OPERATOR_INFO * stiff_info,const OPERATOR_INFO * mass_info,const DOF_REAL_VEC_D * u_h_d)216 fill_sys_info_instat_dow(const OPERATOR_INFO *stiff_info,
217 const OPERATOR_INFO *mass_info,
218 const DOF_REAL_VEC_D *u_h_d)
219 {
220 SYSTEM_FILL_INFO_DOW_INSTAT *info;
221
222 info = MEM_CALLOC(1, SYSTEM_FILL_INFO_DOW_INSTAT);
223
224 info->stiff_info = fill_matrix_info(stiff_info, NULL);
225 info->mass_info = fill_matrix_info(mass_info, NULL);
226
227 info->elsii.krn_blk_type = info->stiff_info->krn_blk_type;
228
229 info->elsii.row_fe_space = info->mass_info->row_fe_space;
230 info->elsii.col_fe_space = info->mass_info->col_fe_space;
231 if (info->elsii.col_fe_space == NULL) {
232 info->elsii.col_fe_space = info->elsii.row_fe_space;
233 }
234
235 info->elsii.el_update_fct = element_system_dow_instat;
236
237 info->n_row = info->elsii.row_fe_space->bas_fcts->n_bas_fcts;
238 info->n_row_max = info->elsii.row_fe_space->bas_fcts->n_bas_fcts_max;
239 info->n_col = info->elsii.col_fe_space->bas_fcts->n_bas_fcts;
240 info->n_col_max = info->elsii.col_fe_space->bas_fcts->n_bas_fcts_max;
241
242 info->elsii.el_matrix = get_el_matrix(info->elsii.row_fe_space,
243 info->elsii.col_fe_space,
244 info->elsii.krn_blk_type);
245 info->elsii.el_load = get_el_real_vec_d(info->elsii.row_fe_space->bas_fcts);
246 info->elsii.u_h_loc = get_el_real_vec_d(info->elsii.col_fe_space->bas_fcts);
247
248 info->elsii.fill_flag =
249 info->stiff_info->fill_flag | info->mass_info->fill_flag;
250 BNDRY_FLAGS_CPY(info->elsii.dirichlet_bndry,
251 info->mass_info->dirichlet_bndry);
252 BNDRY_FLAGS_OR(info->elsii.dirichlet_bndry,
253 info->stiff_info->dirichlet_bndry);
254 if (!BNDRY_FLAGS_IS_INTERIOR(info->elsii.dirichlet_bndry)) {
255 info->elsii.fill_flag |= FILL_BOUND;
256 if (info->elsii.row_fe_space->mesh->is_periodic &&
257 !(info->elsii.row_fe_space->admin->flags & ADM_PERIODIC))
258 info->elsii.fill_flag |= FILL_NON_PERIODIC;
259 }
260
261 info->u_h_d = u_h_d;
262
263 return &info->elsii;
264 }
265
free_sys_info_instat(EL_SYS_INFO_INSTAT * elsii)266 void free_sys_info_instat(EL_SYS_INFO_INSTAT *elsii)
267 {
268 SYSTEM_FILL_INFO_INSTAT *info = (SYSTEM_FILL_INFO_INSTAT *)elsii;
269
270 free_el_matrix(elsii->el_matrix);
271 free_el_real_vec(elsii->el_load);
272 free_el_real_vec((EL_REAL_VEC *)elsii->u_h_loc);
273 MEM_FREE(info->stiff_info, 1, EL_MATRIX_INFO);
274 MEM_FREE(info->mass_info, 1, EL_MATRIX_INFO);
275 MEM_FREE(info, 1, SYSTEM_FILL_INFO_INSTAT);
276 }
277
free_sys_info_dow_instat(EL_SYS_INFO_DOW_INSTAT * elsii)278 void free_sys_info_dow_instat(EL_SYS_INFO_DOW_INSTAT *elsii)
279 {
280 SYSTEM_FILL_INFO_DOW_INSTAT *info = (SYSTEM_FILL_INFO_DOW_INSTAT *)elsii;
281
282 free_el_matrix(elsii->el_matrix);
283 free_el_real_vec_d(elsii->el_load);
284 free_el_real_vec_d((EL_REAL_VEC_D *)elsii->u_h_loc);
285 MEM_FREE(info->stiff_info, 1, EL_MATRIX_INFO);
286 MEM_FREE(info->mass_info, 1, EL_MATRIX_INFO);
287 MEM_FREE(info, 1, SYSTEM_FILL_INFO_DOW_INSTAT);
288 }
289
290 /* Assemble (1/tau M + theta S) and add (1/tau M uh + (1-theta) S uh)
291 * to fh.
292 */
update_system_instat(DOF_MATRIX * dof_matrix,DOF_REAL_VEC * f_h,REAL tau,REAL theta,EL_SYS_INFO_INSTAT * elsii)293 void update_system_instat(DOF_MATRIX *dof_matrix,
294 DOF_REAL_VEC *f_h,
295 REAL tau,
296 REAL theta,
297 EL_SYS_INFO_INSTAT *elsii)
298 {
299 const BAS_FCTS *row_bfcts = elsii->row_fe_space->bas_fcts;
300 const EL_DOF_VEC *row_dof, *col_dof;
301 EL_SCHAR_VEC *bound = NULL;
302 const EL_BNDRY_VEC *bndry_bits;
303 bool use_get_bound;
304
305 BNDRY_FLAGS_CPY(dof_matrix->dirichlet_bndry, elsii->dirichlet_bndry);
306
307 use_get_bound = !BNDRY_FLAGS_IS_INTERIOR(dof_matrix->dirichlet_bndry);
308
309 if (use_get_bound) {
310 bound = get_el_schar_vec(elsii->row_fe_space->bas_fcts);
311 }
312
313 TRAVERSE_FIRST(dof_matrix->row_fe_space->mesh,
314 -1,
315 elsii->fill_flag|CALL_LEAF_EL|FILL_COORDS) {
316
317 if (elsii->el_update_fct(el_info, tau, theta, elsii) == INIT_EL_TAG_NULL) {
318 continue;
319 }
320
321 row_dof = get_dof_indices(NULL, f_h->fe_space, el_info->el);
322 if (elsii->row_fe_space != elsii->col_fe_space) {
323 col_dof =
324 get_dof_indices(NULL, elsii->col_fe_space, el_info->el);
325 } else {
326 col_dof = row_dof;
327 }
328
329 if (use_get_bound) {
330 bndry_bits = get_bound(NULL, row_bfcts, el_info);
331 dirichlet_map(bound, bndry_bits, dof_matrix->dirichlet_bndry);
332 }
333
334 add_element_matrix(dof_matrix, 1.0, elsii->el_matrix, NoTranspose,
335 row_dof, col_dof, bound);
336
337 add_element_vec(f_h, 1.0, elsii->el_load, row_dof, bound);
338
339 } TRAVERSE_NEXT();
340
341 if (use_get_bound) {
342 free_el_schar_vec(bound);
343 }
344 }
345
update_system_instat_dow(DOF_MATRIX * dof_matrix,DOF_REAL_VEC_D * f_h,REAL tau,REAL theta,EL_SYS_INFO_DOW_INSTAT * elsii)346 void update_system_instat_dow(DOF_MATRIX *dof_matrix,
347 DOF_REAL_VEC_D *f_h,
348 REAL tau,
349 REAL theta,
350 EL_SYS_INFO_DOW_INSTAT *elsii)
351 {
352 const BAS_FCTS *row_bfcts = elsii->row_fe_space->bas_fcts;
353 const EL_DOF_VEC *row_dof, *col_dof;
354 EL_SCHAR_VEC *bound = NULL;
355 const EL_BNDRY_VEC *bndry_bits;
356 bool use_get_bound;
357
358 BNDRY_FLAGS_CPY(dof_matrix->dirichlet_bndry, elsii->dirichlet_bndry);
359
360 use_get_bound = !BNDRY_FLAGS_IS_INTERIOR(dof_matrix->dirichlet_bndry);
361
362 if (use_get_bound) {
363 bound = get_el_schar_vec(elsii->row_fe_space->bas_fcts);
364 }
365
366 TRAVERSE_FIRST(dof_matrix->row_fe_space->mesh,
367 -1,
368 elsii->fill_flag|CALL_LEAF_EL|FILL_COORDS) {
369
370 if (elsii->el_update_fct(el_info, tau, theta, elsii) == INIT_EL_TAG_NULL) {
371 continue;
372 }
373
374 row_dof = get_dof_indices(NULL, f_h->fe_space, el_info->el);
375 if (elsii->row_fe_space != elsii->col_fe_space) {
376 col_dof =
377 get_dof_indices(NULL, elsii->col_fe_space, el_info->el);
378 } else {
379 col_dof = row_dof;
380 }
381
382 if (use_get_bound) {
383 bndry_bits = get_bound(NULL, row_bfcts, el_info);
384 dirichlet_map(bound, bndry_bits, dof_matrix->dirichlet_bndry);
385 }
386
387 add_element_matrix(dof_matrix, 1.0, elsii->el_matrix, NoTranspose,
388 row_dof, col_dof, bound);
389
390 add_element_vec_dow(f_h, 1.0, elsii->el_load, row_dof, bound);
391
392 } TRAVERSE_NEXT();
393
394 if (use_get_bound) {
395 free_el_schar_vec(bound);
396 }
397 }
398