1 /*--------------------------------------------------------------------------*/
2 /* ALBERTA: an Adaptive multi Level finite element toolbox using */
3 /* Bisectioning refinement and Error control by Residual */
4 /* Techniques for scientific Applications */
5 /* */
6 /*--------------------------------------------------------------------------*/
7 /* */
8 /* file: nlsolve.c */
9 /* */
10 /* description: solver for an nonlinear elliptic model problem */
11 /* */
12 /* -k \Delta u + \sigma u^4 = f in \Omega */
13 /* u = g on \partial \Omega */
14 /* */
15 /* definitions and settings for different problems */
16 /* */
17 /*--------------------------------------------------------------------------*/
18 /* */
19 /* authors: Alfred Schmidt */
20 /* Zentrum fuer Technomathematik */
21 /* Fachbereich 3 Mathematik/Informatik */
22 /* Universitaet Bremen */
23 /* Bibliothekstr. 2 */
24 /* D-28359 Bremen, Germany */
25 /* */
26 /* Kunibert G. Siebert */
27 /* Institut fuer Mathematik */
28 /* Universitaet Augsburg */
29 /* Universitaetsstr. 14 */
30 /* D-86159 Augsburg, Germany */
31 /* */
32 /* http://www.alberta-fem.de/ */
33 /* */
34 /*--------------------------------------------------------------------------*/
35 /* (c) by A. Schmidt and K.G. Siebert (1996-2004) */
36 /* */
37 /* This program is free software; you can redistribute it and/or modify */
38 /* it under the terms of the GNU General Public License as published by */
39 /* the Free Software Foundation; either version 2 of the License, or */
40 /* any later version. */
41 /* */
42 /* This program is distributed in the hope that it will be useful, */
43 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
44 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
45 /* GNU General Public License for more details. */
46 /*--------------------------------------------------------------------------*/
47
48 #include "nonlin.h"
49
50 static REAL k = 1.0, sigma = 1.0;
51
52 /*--------------------------------------------------------------------------*/
53 /* Problem 0: problem with exact solution */
54 /* solution is the `exponential peak' */
55 /*--------------------------------------------------------------------------*/
56
u_0(const REAL_D x)57 static REAL u_0(const REAL_D x)
58 {
59 REAL x2 = SCP_DOW(x,x);
60 return(exp(-10.0*x2));
61 }
62
grd_u_0(const REAL_D x)63 static const REAL *grd_u_0(const REAL_D x)
64 {
65 static REAL_D grd;
66 REAL ux = exp(-10.0*SCP_DOW(x,x));
67 int n;
68
69 for (n = 0; n < DIM_OF_WORLD; n++)
70 grd[n] = -20.0*x[n]*ux;
71
72 return(grd);
73 }
74
f_0(const REAL_D x)75 static REAL f_0(const REAL_D x)
76 {
77 REAL r2 = SCP_DOW(x,x), ux = exp(-10.0*r2), ux4 = ux*ux*ux*ux;
78 return(sigma*ux4 - k*(400.0*r2 - 20.0*DIM)*ux);
79 }
80
81 /*--------------------------------------------------------------------------*/
82 /* Problem 1: */
83 /* compute a stable and an unstable solution, depending on the inital */
84 /* choice of u_h */
85 /* there is no true solution */
86 /*--------------------------------------------------------------------------*/
87
88 static REAL U0 = 0.0;
89
g_1(const REAL_D x)90 static REAL g_1(const REAL_D x)
91 {
92 #if DIM_OF_WORLD == 1
93 return(4.0*U0*x[0]*(1.0-x[0]));
94 #endif
95 #if DIM_OF_WORLD == 2
96 return(16.0*U0*x[0]*(1.0-x[0])*x[1]*(1.0-x[1]));
97 #endif
98 #if DIM_OF_WORLD == 3
99 return(64.0*U0*x[0]*(1.0-x[0])*x[1]*(1.0-x[1])*x[2]*(1.0-x[2]));
100 #endif
101 }
102
f_1(const REAL_D x)103 static REAL f_1(const REAL_D x)
104 {
105 return(1.0);
106 }
107
108 /*--------------------------------------------------------------------------*/
109 /* Problem 2: */
110 /* physical problem: heat transport in a thin plate */
111 /* no true solution known */
112 /*--------------------------------------------------------------------------*/
113
114
115 static REAL g2 = 300.0, sigma_uext4 = 0.0;
g_2(const REAL_D x)116 static REAL g_2(const REAL_D x)
117 {
118 return(g2);
119 }
120
f_2(const REAL_D x)121 static REAL f_2(const REAL_D x)
122 {
123 if (x[0] >= -0.25 && x[0] <= 0.25 && x[1] >= -0.25 && x[1] <= 0.25)
124 return(150.0 + sigma_uext4);
125 else
126 return(sigma_uext4);
127 }
128
129 /*--------------------------------------------------------------------------*/
130 /* init_problem: */
131 /* initialize k and sigma, adjust function pointers and the read macro */
132 /* triangulation */
133 /*--------------------------------------------------------------------------*/
134
init_problem(MESH * mesh)135 const PROB_DATA *init_problem(MESH *mesh)
136 {
137 FUNCNAME("init_problem");
138 static PROB_DATA prob_data = {0};
139 int pn = 2, n_refine = 0;
140
141 GET_PARAMETER(1, "problem number", "%d", &pn);
142 switch (pn)
143 {
144 case 0: /*--- problem with known true solution -----------------------*/
145 k = 1.0;
146 sigma = 1.0;
147
148 prob_data.g = u_0;
149 prob_data.f = f_0;
150
151 prob_data.u = u_0;
152 prob_data.grd_u = grd_u_0;
153
154 read_macro(mesh, "Macro/macro-big.amc", nil);
155 break;
156 case 1: /*--- problem for computing a stable and an unstable sol. ----*/
157 k = 1.0;
158 sigma = 1.0;
159
160 prob_data.g = g_1;
161 prob_data.f = f_1;
162
163 prob_data.u0 = g_1;
164 GET_PARAMETER(1, "U0", "%f", &U0);
165
166 read_macro(mesh, "Macro/macro.amc", nil);
167 break;
168 case 2: /*--- physical problem ---------------------------------------*/
169 TEST_EXIT(DIM==2)("problem 2 makes sense only in 2d!\n");
170 k = 2.0;
171 sigma = 5.67e-8;
172 sigma_uext4 = sigma*273*273*273*273;
173
174 prob_data.g = g_2;
175 prob_data.f = f_2;
176 read_macro(mesh, "Macro/macro-big.amc", nil);
177 break;
178 default:
179 ERROR_EXIT("no problem defined with problem no. %d\n", pn);
180 }
181 prob_data.k = k;
182 prob_data.sigma = sigma;
183
184 GET_PARAMETER(1, "global refinements", "%d", &n_refine);
185 global_refine(mesh, n_refine*DIM);
186
187 return(&prob_data);
188 }
189