1 /*
2 * ***************************************************************************
3 * MC = < Manifold Code >
4 * Copyright (C) 1994-- Michael Holst
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 * rcsid="$Id: pde.c,v 1.8 2010/08/12 05:19:16 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: pde.c
27 *
28 * Purpose: Class PDE: methods.
29 *
30 * Author: Michael Holst
31 * ***************************************************************************
32 */
33
34 #include "pde_p.h"
35
36 VEMBED(rcsid="$Id: pde.c,v 1.8 2010/08/12 05:19:16 fetk Exp $")
37
38 /*
39 * ***************************************************************************
40 * Class PDE: Inlineable methods
41 * ***************************************************************************
42 */
43 #if !defined(VINLINE_PDE)
44
45 /*
46 * ***************************************************************************
47 * Routine: PDE_setDim
48 *
49 * Purpose: Set the extrinsic spatial dimension.
50 *
51 * Author: Michael Holst
52 * ***************************************************************************
53 */
PDE_setDim(PDE * thee,int d)54 VPUBLIC void PDE_setDim(PDE *thee, int d)
55 {
56 thee->dim = d;
57 }
58
59 /*
60 * ***************************************************************************
61 * Routine: PDE_setDimII
62 *
63 * Purpose: Set the intrinsic spatial dimension.
64 *
65 * Author: Michael Holst
66 * ***************************************************************************
67 */
PDE_setDimII(PDE * thee,int d)68 VPUBLIC void PDE_setDimII(PDE *thee, int d)
69 {
70 thee->dimII = d;
71 }
72
73 /*
74 * ***************************************************************************
75 * Routine: PDE_dim
76 *
77 * Purpose: Return the extrinsic spatial dimension.
78 *
79 * Author: Michael Holst
80 * ***************************************************************************
81 */
PDE_dim(PDE * thee)82 VPUBLIC int PDE_dim(PDE *thee)
83 {
84 return thee->dim;
85 }
86
87 /*
88 * ***************************************************************************
89 * Routine: PDE_dimII
90 *
91 * Purpose: Return the intrinsic spatial dimension.
92 *
93 * Author: Michael Holst
94 * ***************************************************************************
95 */
PDE_dimII(PDE * thee)96 VPUBLIC int PDE_dimII(PDE *thee)
97 {
98 return thee->dimII;
99 }
100
101 /*
102 * ***************************************************************************
103 * Routine: PDE_vec
104 *
105 * Purpose: Return the PDE product space dimension.
106 *
107 * Author: Michael Holst
108 * ***************************************************************************
109 */
PDE_vec(PDE * thee)110 VPUBLIC int PDE_vec(PDE *thee)
111 {
112 return thee->vec;
113 }
114
115 #endif /* if !defined(VINLINE_PDE) */
116 /*
117 * ***************************************************************************
118 * Class PDE: Non-inlineable methods
119 * ***************************************************************************
120 */
121
122 /*
123 * ***************************************************************************
124 * Routine: bisectEdge_default
125 *
126 * Purpose: DEFAULT: Define the way manifold edges are bisected.
127 *
128 * Author: Michael Holst
129 * ***************************************************************************
130 */
bisectEdge_default(int dim,int dimII,int edgeType,int chart[],double vx[][3])131 VPRIVATE void bisectEdge_default(int dim, int dimII,
132 int edgeType, int chart[], double vx[][3])
133 {
134 int i;
135 for (i=0; i<dimII; i++)
136 vx[2][i] = .5 * (vx[0][i] + vx[1][i]);
137 }
138
139 /*
140 * ***************************************************************************
141 * Routine: mapBoundary_default
142 *
143 * Purpose: DEFAULT: Define the way the manifold boundary is treated.
144 *
145 * Author: Michael Holst
146 * ***************************************************************************
147 */
mapBoundary_default(int dim,int dimII,int vertexType,int chart,double vx[3])148 VPRIVATE void mapBoundary_default(int dim, int dimII,
149 int vertexType, int chart, double vx[3])
150 {
151 }
152
153 /*
154 * ***************************************************************************
155 * Routine: markSimplex_default
156 *
157 * Purpose: DEFAULT: Simplex marking routine for refinement.
158 *
159 * Author: Michael Holst
160 * ***************************************************************************
161 */
markSimplex_default(int dim,int dimII,int simplexType,int faceType[4],int vertexType[4],int chart[],double vx[][3],void * data)162 VPRIVATE int markSimplex_default(int dim, int dimII,
163 int simplexType, int faceType[4], int vertexType[4],
164 int chart[], double vx[][3], void *data)
165 {
166 int j, k, less, more;
167 double radius, d[4];
168
169 /* radius = radius of a refinement sphere for testing */
170 radius = 0.1; /* must be > 0 */
171 less = 0;
172 more = 0;
173 for (j=0; j<dim+1; j++) {
174 d[j] = 0.0;
175 for (k=0; k<3; k++)
176 d[j] += VSQR( vx[j][k] );
177 d[j] = VSQRT( d[j] );
178 if (d[j] <= radius+VSMALL) less = 1;
179 else more = 1;
180 }
181
182 /* return true if simplex touches or stradles surface of sphere */
183 return ( less && more );
184 }
185
186 /*
187 * ***************************************************************************
188 * Routine: oneChart_default
189 *
190 * Purpose: DEFAULT: Select a single unified chart for a set of two or
191 * more vertices whose coordinates may be given with respect to
192 * different charts. Then transform all of the coordinates of
193 * the vertex set to be with respect to the single selected
194 * "unified" chart.
195 *
196 * Author: Michael Holst
197 * ***************************************************************************
198 */
oneChart_default(int dim,int dimII,int objType,int chart[],double vx[][3],int dimV)199 VPRIVATE void oneChart_default(int dim, int dimII, int objType,
200 int chart[], double vx[][3], int dimV)
201 {
202 VASSERT( (2 <= dim) && (dim <= 3) );
203 VASSERT( (2 <= dimII) && (dimII <= 3) );
204 VASSERT( (1 <= dimV) && (dimV <= 4) );
205 VASSERT( (0 <= objType) );
206 }
207
208 /*
209 * ***************************************************************************
210 * Routine: PDE_ctor_default
211 *
212 * Purpose: Construct a fake differential equation object in case there
213 * is not one provided.
214 *
215 * Author: Michael Holst
216 * ***************************************************************************
217 */
PDE_ctor_default(void)218 VPUBLIC PDE* PDE_ctor_default(void)
219 {
220 int i;
221 PDE *thee = VNULL;
222
223 /* create some space for the pde object */
224 thee = Vmem_malloc( VNULL, 1, sizeof(PDE) );
225
226 /* PDE-specific parameters and function pointers */
227 thee->initAssemble = VNULL; /* once-per-assembly initialization */
228 thee->initElement = VNULL; /* once-per-element initialization */
229 thee->initFace = VNULL; /* once-per-face initialization */
230 thee->initPoint = VNULL; /* once-per-point initialization */
231 thee->Fu = VNULL; /* nonlinear strong form */
232 thee->Ju = VNULL; /* nonlinear energy functional */
233 thee->Fu_v = VNULL; /* nonlinear weak form */
234 thee->DFu_wv = VNULL; /* bilinear linearization weak form */
235 thee->p_wv = VNULL; /* bilinear mass density form */
236 thee->delta = VNULL; /* delta function source term */
237 thee->u_D = VNULL; /* dirichlet func and initial guess */
238 thee->u_T = VNULL; /* analytical soln for testing */
239 thee->vec = 1; /* unknowns per spatial point; */
240 thee->sym[0][0] = 0; /* symmetries of bilinear form */
241 thee->est[0] = 1.0; /* error estimator weights */
242 for (i=0; i<VMAX_BDTYPE; i++) /* boundary type remappings */
243 thee->bmap[0][i] = i;
244
245 /* Manifold-specific function pointers */
246 thee->bisectEdge = bisectEdge_default; /* edge bisection rule */
247 thee->mapBoundary = mapBoundary_default; /* boundary recovery rule */
248 thee->markSimplex = markSimplex_default; /* simplex marking rule */
249 thee->oneChart = oneChart_default; /* coordinate transformations */
250
251 /* Element-specific function pointers */
252 thee->simplexBasisInit = VNULL; /* initialization of bases */
253 thee->simplexBasisForm = VNULL; /* form trial & test bases */
254
255 /* return the new pde object */
256 return thee;
257 }
258
259 /*
260 * ***************************************************************************
261 * Routine: PDE_dtor_default
262 *
263 * Purpose: Destroy a fake differential equation object.
264 *
265 * Author: Michael Holst
266 * ***************************************************************************
267 */
PDE_dtor_default(PDE ** thee)268 VPUBLIC void PDE_dtor_default(PDE **thee)
269 {
270 VASSERT( (*thee) != VNULL );
271 if ((*thee) != VNULL) {
272
273 Vmem_free( VNULL, 1, sizeof(PDE), (void**)thee );
274
275 (*thee) = VNULL;
276 }
277 }
278
279 /*
280 * ***************************************************************************
281 * Routine: PDE_sym
282 *
283 * Purpose: Return the bilinear form minor symmetry (i,j).
284 *
285 * Author: Michael Holst
286 * ***************************************************************************
287 */
PDE_sym(PDE * thee,int i,int j)288 VPUBLIC int PDE_sym(PDE *thee, int i, int j)
289 {
290 return thee->sym[i][j];
291 }
292
293