183{
184
185
186
187
188
189
190
191 if (type == "Bifurcation")
192 {
193 NekMatrix<NekDouble> J(6, 6);
194 Array<OneD, NekDouble> c(3);
195
196 for (int i = 0; i < 3; ++i)
197 {
198 GetC(c[i],
beta[i], Au[i], A0[i], alpha[i]);
199 }
200
201 J.SetValue(0, 0, 1);
202 J.SetValue(0, 1, 0);
203 J.SetValue(0, 2, 0);
204 J.SetValue(0, 3, c[0] / Au[0]);
205 J.SetValue(0, 4, 0);
206 J.SetValue(0, 5, 0);
207
208 J.SetValue(1, 0, 0);
209 J.SetValue(1, 1, 1);
210 J.SetValue(1, 2, 0);
211 J.SetValue(1, 3, 0);
212 J.SetValue(1, 4, -c[1] / Au[1]);
213 J.SetValue(1, 5, 0);
214
215 J.SetValue(2, 0, 0);
216 J.SetValue(2, 1, 0);
217 J.SetValue(2, 2, 1);
218 J.SetValue(2, 3, 0);
219 J.SetValue(2, 4, 0);
220 J.SetValue(2, 5, -c[2] / Au[2]);
221
222 J.SetValue(3, 0, Au[0]);
223 J.SetValue(3, 1, -Au[1]);
224 J.SetValue(3, 2, -Au[2]);
225 J.SetValue(3, 3, uu[0]);
226 J.SetValue(3, 4, -uu[1]);
227 J.SetValue(3, 5, -uu[2]);
228
229 J.SetValue(4, 0, 2 * uu[0]);
230 J.SetValue(4, 1, -2 * uu[1]);
231 J.SetValue(4, 2, 0);
232 J.SetValue(4, 3, 2 * c[0] * c[0] / Au[0]);
233 J.SetValue(4, 4, -2 * c[1] * c[1] / Au[1]);
234 J.SetValue(4, 5, 0);
235
236 J.SetValue(5, 0, 2 * uu[0]);
237 J.SetValue(5, 1, 0);
238 J.SetValue(5, 2, -2 * uu[2]);
239 J.SetValue(5, 3, 2 * c[0] * c[0] / Au[0]);
240 J.SetValue(5, 4, 0);
241 J.SetValue(5, 5, -2 * c[2] * c[2] / Au[2]);
242
243 invJ = J;
244 invJ.Invert();
245 }
246 else if (type == "Merge")
247 {
248 NekMatrix<NekDouble> J(6, 6);
249 Array<OneD, NekDouble> c(3);
250
251 for (int i = 0; i < 3; ++i)
252 {
253 GetC(c[i],
beta[i], Au[i], A0[i], alpha[i]);
254 }
255
256 J.SetValue(0, 0, 1);
257 J.SetValue(0, 1, 0);
258 J.SetValue(0, 2, 0);
259 J.SetValue(0, 3, -c[0] / Au[0]);
260 J.SetValue(0, 4, 0);
261 J.SetValue(0, 5, 0);
262
263 J.SetValue(1, 0, 0);
264 J.SetValue(1, 1, 1);
265 J.SetValue(1, 2, 0);
266 J.SetValue(1, 3, 0);
267 J.SetValue(1, 4, c[1] / Au[1]);
268 J.SetValue(1, 5, 0);
269
270 J.SetValue(2, 0, 0);
271 J.SetValue(2, 1, 0);
272 J.SetValue(2, 2, 1);
273 J.SetValue(2, 3, 0);
274 J.SetValue(2, 4, 0);
275 J.SetValue(2, 5, c[2] / Au[2]);
276
277 J.SetValue(3, 0, Au[0]);
278 J.SetValue(3, 1, -Au[1]);
279 J.SetValue(3, 2, -Au[2]);
280 J.SetValue(3, 3, uu[0]);
281 J.SetValue(3, 4, -uu[1]);
282 J.SetValue(3, 5, -uu[2]);
283
284 J.SetValue(4, 0, 2 * uu[0]);
285 J.SetValue(4, 1, -2 * uu[1]);
286 J.SetValue(4, 2, 0);
287 J.SetValue(4, 3, 2 * c[0] * c[0] / Au[0]);
288 J.SetValue(4, 4, -2 * c[1] * c[1] / Au[1]);
289 J.SetValue(4, 5, 0);
290
291 J.SetValue(5, 0, 2 * uu[0]);
292 J.SetValue(5, 1, 0);
293 J.SetValue(5, 2, -2 * uu[2]);
294 J.SetValue(5, 3, 2 * c[0] * c[0] / Au[0]);
295 J.SetValue(5, 4, 0);
296 J.SetValue(5, 5, -2 * c[2] * c[2] / Au[2]);
297
298 invJ = J;
299 invJ.Invert();
300 }
301 else if (type == "Interface")
302 {
303 NekMatrix<NekDouble> J(4, 4);
304 Array<OneD, NekDouble> c(2);
305
306 for (int i = 0; i < 2; ++i)
307 {
308 GetC(c[i],
beta[i], Au[i], A0[i], alpha[i]);
309 }
310
311 J.SetValue(0, 0, 1);
312 J.SetValue(0, 1, 0);
313 J.SetValue(0, 2, c[0] / Au[0]);
314 J.SetValue(0, 3, 0);
315
316 J.SetValue(1, 0, 0);
317 J.SetValue(1, 1, 1);
318 J.SetValue(1, 2, 0);
319 J.SetValue(1, 3, -c[1] / Au[1]);
320
321 J.SetValue(2, 0, Au[0]);
322 J.SetValue(2, 1, -Au[1]);
323 J.SetValue(2, 2, uu[0]);
324 J.SetValue(2, 3, -uu[1]);
325
326 J.SetValue(3, 0, 2 * uu[0]);
327 J.SetValue(3, 1, -2 * uu[1]);
328 J.SetValue(3, 2, 2 * c[0] * c[0] / Au[0]);
329 J.SetValue(3, 3, -2 * c[1] * c[1] / Au[1]);
330
331 invJ = J;
332 invJ.Invert();
333 }
334}