Evaluator.hxx 36.2 KB
Newer Older
Michael Turner's avatar
Michael Turner committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
////////////////////////////////////////////////////////////////////////////////
//
//  File: ProcessJac.h
//
//  For more information, please see: http://www.nektar.info/
//
//  The MIT License
//
//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
//  Department of Aeronautics, Imperial College London (UK), and Scientific
//  Computing and Imaging Institute, University of Utah (USA).
//
//  License for the specific language governing rights and limitations under
//  Permission is hereby granted, free of charge, to any person obtaining a
//  copy of this software and associated documentation files (the "Software"),
//  to deal in the Software without restriction, including without limitation
//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
//  and/or sell copies of the Software, and to permit persons to whom the
//  Software is furnished to do so, subject to the following conditions:
//
//  The above copyright notice and this permission notice shall be included
//  in all copies or substantial portions of the Software.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
//  DEALINGS IN THE SOFTWARE.
//
//  Description: Calculate jacobians of elements.
//
////////////////////////////////////////////////////////////////////////////////

#ifndef UTILITIES_NEKMESH_NODEOPTI_EVALUATOR
#define UTILITIES_NEKMESH_NODEOPTI_EVALUATOR

namespace Nektar
{
namespace Utilities
{

Michael Turner's avatar
Michael Turner committed
44
45
using namespace std;

Michael Turner's avatar
Michael Turner committed
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
template<int DIM> inline NekDouble Determinant(NekDouble jac[DIM][DIM])
{
    return 0.0;
}

template<> inline NekDouble Determinant<2>(NekDouble jac[2][2])
{
    return jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0];
}

template<> inline NekDouble Determinant<3>(NekDouble jac[3][3])
{
    return jac[0][0] * (jac[1][1]*jac[2][2] - jac[2][1]*jac[1][2])
          -jac[0][1] * (jac[1][0]*jac[2][2] - jac[1][2]*jac[2][0])
          +jac[0][2] * (jac[1][0]*jac[2][1] - jac[1][1]*jac[2][0]);
}

template<int DIM> inline void InvTrans(NekDouble in[DIM][DIM],
                                       NekDouble out[DIM][DIM])
{
}

template<>
inline void InvTrans<2>(NekDouble in[2][2], NekDouble out[2][2])
{
    NekDouble invDet = 1.0 / Determinant(in);
    out[0][0] =  in[1][1] * invDet;
    out[1][0] = -in[0][1] * invDet;
    out[0][1] = -in[1][0] * invDet;
    out[1][1] =  in[0][0] * invDet;
}

template<>
inline void InvTrans<3>(NekDouble in[3][3], NekDouble out[3][3])
{
    NekDouble invdet = 1.0 / Determinant(in);
    out[0][0] =  (in[1][1]*in[2][2]-in[2][1]*in[1][2])*invdet;
    out[1][0] = -(in[0][1]*in[2][2]-in[0][2]*in[2][1])*invdet;
    out[2][0] =  (in[0][1]*in[1][2]-in[0][2]*in[1][1])*invdet;
    out[0][1] = -(in[1][0]*in[2][2]-in[1][2]*in[2][0])*invdet;
    out[1][1] =  (in[0][0]*in[2][2]-in[0][2]*in[2][0])*invdet;
    out[2][1] = -(in[0][0]*in[1][2]-in[1][0]*in[0][2])*invdet;
    out[0][2] =  (in[1][0]*in[2][1]-in[2][0]*in[1][1])*invdet;
    out[1][2] = -(in[0][0]*in[2][1]-in[2][0]*in[0][1])*invdet;
    out[2][2] =  (in[0][0]*in[1][1]-in[1][0]*in[0][1])*invdet;
}

93
94
95
96
97
98
99
100
template<int DIM> inline void EMatrix(NekDouble in[DIM][DIM],
                                       NekDouble out[DIM][DIM])
{
}

template<>
inline void EMatrix<2>(NekDouble in[2][2], NekDouble out[2][2])
{
101
102
103
104
    out[0][0] =  0.5*(in[0][0]*in[0][0] + in[1][0]*in[1][0] - 1.0);
    out[1][0] =  0.5*(in[0][0]*in[0][1] + in[1][0]*in[1][1]);
    out[0][1] =  0.5*(in[0][0]*in[0][1] + in[1][0]*in[1][1]);
    out[1][1] =  0.5*(in[1][1]*in[1][1] + in[0][1]*in[0][1] - 1.0);
105
106
107
108
109
}

template<>
inline void EMatrix<3>(NekDouble in[3][3], NekDouble out[3][3])
{
110
111
    out[0][0] =  0.5*(in[0][0]*in[0][0] + in[1][0]*in[1][0] + in[2][0]*in[2][0]- 1.0);
    out[1][0] =  0.5*(in[0][0]*in[1][0] + in[1][0]*in[1][1] + in[2][0]*in[2][1]);
112
    out[0][1] =  out[1][0];
113
    out[2][0] =  0.5*(in[0][0]*in[0][2] + in[1][0]*in[1][2] + in[2][0]*in[2][2]);
114
    out[0][2] =  out[2][0];
115
116
    out[1][1] =  0.5*(in[0][1]*in[0][1] + in[1][1]*in[1][1] + in[2][1]*in[2][1]- 1.0);
    out[1][2] =  0.5*(in[0][1]*in[0][2] + in[1][1]*in[1][2] + in[2][1]*in[2][2]);
117
    out[2][1] =  out[1][2];
118
119
120
    out[2][2] =  0.5*(in[0][2]*in[0][2] + in[1][2]*in[1][2] + in[2][2]*in[2][2]- 1.0);
}

121
122
123
template<int DIM> inline void LEM2(NekDouble jacIdeal[DIM][DIM],
                                   NekDouble jacDerivPhi[DIM][DIM][DIM],
                                   NekDouble ret[DIM][DIM][DIM])
Michael Turner's avatar
Michael Turner committed
124
{
125
    for(int k = 0; k < DIM; k++)
126
    {
127
128
        NekDouble part1[DIM][DIM], part2[DIM][DIM];
        for (int m = 0; m < DIM; ++m)
129
        {
130
            for (int n = 0; n < DIM; ++n)
131
            {
132
133
134
135
136
137
138
                part1[m][n] = 0.0;
                part2[m][n] = 0.0;
                for (int l = 0; l < DIM; ++l)
                {
                    part1[m][n] += jacDerivPhi[k][l][m] * jacIdeal[l][n];
                    part2[m][n] += jacIdeal[l][m] * jacDerivPhi[k][l][n];
                }
139
140
141
            }
        }

142
        for (int m = 0; m < DIM; ++m)
143
        {
144
            for (int n = 0; n < DIM; ++n)
145
            {
146
                ret[k][m][n] = 0.5*(part1[m][n] + part2[m][n]);
147
148
149
            }
        }
    }
150
}
151

152
153
154
155
template<int DIM> inline void LEM3(NekDouble jacDerivPhi[DIM][DIM][DIM],
                                   NekDouble ret[DIM][DIM][DIM][DIM])
{
    for(int j = 0; j < DIM; j++)
156
    {
157
        for(int k = 0; k < DIM; k++)
158
        {
159
160
            NekDouble part1[DIM][DIM], part2[DIM][DIM];
            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
161
            {
162
163
164
165
166
167
                for (int n = 0; n < DIM; ++n)
                {
                    part1[m][n] = 0.0;
                    part2[m][n] = 0.0;
                    for (int l = 0; l < DIM; ++l)
                    {
168
169
                        part1[m][n] += jacDerivPhi[j][l][m] * jacDerivPhi[k][l][n];
                        part2[m][n] += jacDerivPhi[k][l][m] * jacDerivPhi[j][l][n];
170
171
                    }
                }
Michael Turner's avatar
Michael Turner committed
172
            }
173
174

            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
175
            {
176
177
                for (int n = 0; n < DIM; ++n)
                {
178
                    ret[j][k][m][n] = 0.5*(part1[m][n] + part2[m][n]);
179
                }
Michael Turner's avatar
Michael Turner committed
180
            }
181
182
183
184
        }
    }
}

Michael Turner's avatar
Michael Turner committed
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
template<int DIM> inline NekDouble FrobProd(NekDouble in1[DIM][DIM],
                                            NekDouble in2[DIM][DIM])
{
    NekDouble ret = 0;
    for (int n = 0; n < DIM; ++n)
    {
        for (int l = 0; l < DIM; ++l)
        {
            ret += in1[n][l] * in2[n][l];
        }
    }
    return ret;
}

// Typedef for derivative storage, we use boost::multi_array so we can pass this
// to functions easily
typedef boost::multi_array<NekDouble, 4> DerivArray;

template<int DIM>
inline NekDouble CalcIdealJac(int elmt,
                              int point,
                              DerivArray &deriv,
                              std::vector<ElUtilSharedPtr> &data,
                              NekDouble jacIdeal[DIM][DIM])
{
    for (int m = 0; m < DIM; ++m)
    {
        for (int n = 0; n < DIM; ++n)
        {
            jacIdeal[n][m] = 0.0;
            for (int l = 0; l < DIM; ++l)
            {
                jacIdeal[n][m] += deriv[l][elmt][n][point] *
                    data[elmt]->maps[point][m * 3 + l];
            }
        }
    }

    return Determinant(jacIdeal);
}

template<int DIM>
inline NekDouble FrobeniusNorm(NekDouble inarray[DIM][DIM])
{
    NekDouble ret = 0.0, *start = &inarray[0][0];
    for (int i = 0; i < DIM*DIM; ++i, ++start)
    {
        ret += (*start) * (*start);
    }
    return ret;
}

template<int DIM>
Michael Turner's avatar
Michael Turner committed
238
NekDouble NodeOpti::GetFunctional(NekDouble &minJacNew, bool gradient)
Michael Turner's avatar
Michael Turner committed
239
{
Michael Turner's avatar
Michael Turner committed
240
    map<LibUtilities::ShapeType,vector<ElUtilSharedPtr> >::iterator typeIt;
Michael Turner's avatar
Michael Turner committed
241
    map<LibUtilities::ShapeType,DerivArray> derivs;
Michael Turner's avatar
Michael Turner committed
242

Michael Turner's avatar
Michael Turner committed
243
    for(typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
Michael Turner's avatar
Michael Turner committed
244
    {
245
        const int nElmt  = typeIt->second.size();
Michael Turner's avatar
Michael Turner committed
246
        const int totpts = m_derivUtils[typeIt->first]->ptsStd * nElmt;
247
248
249
250
        NekDouble X[DIM * totpts];

        // Store x/y components of each element sequentially in memory
        for (int i = 0, cnt = 0; i < nElmt; ++i)
Michael Turner's avatar
Michael Turner committed
251
        {
Michael Turner's avatar
Michael Turner committed
252
            for (int j = 0; j < m_derivUtils[typeIt->first]->ptsStd; ++j)
Michael Turner's avatar
Michael Turner committed
253
            {
254
255
                for (int d = 0; d < DIM; ++d)
                {
Michael Turner's avatar
Michael Turner committed
256
                    X[cnt + d*m_derivUtils[typeIt->first]->ptsStd + j] =
257
258
                                        *(typeIt->second[i]->nodes[j][d]);
                }
Michael Turner's avatar
Michael Turner committed
259
            }
Michael Turner's avatar
Michael Turner committed
260
            cnt += DIM*m_derivUtils[typeIt->first]->ptsStd;
Michael Turner's avatar
Michael Turner committed
261
262
        }

263
264
265
266
267
        // Storage for derivatives, ordered by:
        //   - standard coordinate direction
        //   - number of elements
        //   - cartesian coordinate direction
        //   - quadrature points
Michael Turner's avatar
Michael Turner committed
268
        //DerivArray deriv(boost::extents[DIM][nElmt][DIM][derivUtil[typeIt->first]->ptsHigh]);
269
        derivs.insert(std::make_pair(
270
                          typeIt->first,
Michael Turner's avatar
Michael Turner committed
271
                          DerivArray(boost::extents[DIM][nElmt][DIM][m_derivUtils[typeIt->first]->pts])));
Michael Turner's avatar
Michael Turner committed
272

273
274
275
276
        // Calculate x- and y-gradients
        for (int d = 0; d < DIM; ++d)
        {
            Blas::Dgemm(
Michael Turner's avatar
Michael Turner committed
277
278
279
280
281
282
283
284
                'N', 'N',
                m_derivUtils[typeIt->first]->pts, DIM * nElmt,
                m_derivUtils[typeIt->first]->ptsStd, 1.0,
                m_derivUtils[typeIt->first]->VdmD[d].GetRawPtr(),
                m_derivUtils[typeIt->first]->pts, X,
                m_derivUtils[typeIt->first]->ptsStd, 0.0,
                &derivs[typeIt->first][d][0][0][0],
                m_derivUtils[typeIt->first]->pts);
285
        }
Michael Turner's avatar
Michael Turner committed
286
287
    }

288
    minJacNew = std::numeric_limits<double>::max();
Michael Turner's avatar
Michael Turner committed
289
    NekDouble integral = 0.0;
Michael Turner's avatar
Michael Turner committed
290
291
    NekDouble ep = m_minJac < 0.0 ? sqrt(1e-8 + 0.04*m_minJac*m_minJac) : 1e-4;
    //NekDouble ep = m_minJac < 0.0 ? sqrt(gam*(gam-m_minJac)) : gam;
292
    //NekDouble ep = minJac < 0.0 ? sqrt(gam*gam + minJac*minJac) : gam;
Michael Turner's avatar
Michael Turner committed
293
    //NekDouble ep = 1e-2;
Michael Turner's avatar
Michael Turner committed
294
    NekDouble jacIdeal[DIM][DIM], jacDet;
Michael Turner's avatar
Michael Turner committed
295
    m_grad = Array<OneD, NekDouble>(DIM == 2 ? 5 : 9, 0.0);
Michael Turner's avatar
Michael Turner committed
296

Michael Turner's avatar
Michael Turner committed
297
    switch(m_opti)
Michael Turner's avatar
Michael Turner committed
298
299
300
    {
        case eLinEl:
        {
301
            const NekDouble nu = 0.45;
Michael Turner's avatar
Michael Turner committed
302
303
304
            const NekDouble mu = 1.0 / 2.0 / (1.0+nu);
            const NekDouble K  = 1.0 / 3.0 / (1.0 - 2.0 * nu);

Michael Turner's avatar
Michael Turner committed
305
            for(typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
Michael Turner's avatar
Michael Turner committed
306
            {
Michael Turner's avatar
Michael Turner committed
307
                NekVector<NekDouble> &quadW = m_derivUtils[typeIt->first]->quadW;
Michael Turner's avatar
Michael Turner committed
308
                for (int i = 0; i < typeIt->second.size(); ++i)
Michael Turner's avatar
Michael Turner committed
309
                {
Michael Turner's avatar
Michael Turner committed
310
                    for(int k = 0; k < m_derivUtils[typeIt->first]->pts; ++k)
311
                    {
Michael Turner's avatar
Michael Turner committed
312
                        jacDet = CalcIdealJac(i, k, derivs[typeIt->first], typeIt->second, jacIdeal);
313
                        minJacNew = min(minJacNew,jacDet);
314

315
316
317
318
                        NekDouble Emat[DIM][DIM];
                        EMatrix<DIM>(jacIdeal,Emat);

                        NekDouble trEtE = FrobProd<DIM>(Emat,Emat);
319
320
                        NekDouble sigma =
                            0.5*(jacDet + sqrt(jacDet*jacDet + 4.0*ep*ep));
Michael Turner's avatar
Michael Turner committed
321

322
323
324
325
                        if(sigma < numeric_limits<double>::min() && !gradient)
                        {
                            return numeric_limits<double>::max();
                        }
Michael Turner's avatar
Michael Turner committed
326
327
328
329
330
                        ASSERTL0(sigma > numeric_limits<double>::min(),
                                 std::string("dividing by zero ") +
                                 boost::lexical_cast<string>(sigma) + " " +
                                 boost::lexical_cast<string>(jacDet) + " " +
                                 boost::lexical_cast<string>(ep));
Michael Turner's avatar
Michael Turner committed
331

332
                        NekDouble lsigma = log(sigma);
Michael Turner's avatar
Michael Turner committed
333
                        integral += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) *
334
                                    (K * 0.5 * lsigma * lsigma + mu * trEtE);
Michael Turner's avatar
Michael Turner committed
335

336
                        if(gradient)
337
                        {
338
                            NekDouble jacInvTrans[DIM][DIM];
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
                            NekDouble jacDetDeriv[DIM];

                            NekDouble phiM[DIM][DIM];
                            for (int m = 0; m < DIM; ++m)
                            {
                                for (int n = 0; n < DIM; ++n)
                                {
                                    phiM[n][m] = derivs[typeIt->first][m][i][n][k];
                                }
                            }

                            InvTrans<DIM>(phiM, jacInvTrans);
                            NekDouble derivDet = Determinant<DIM>(phiM);

                            NekDouble basisDeriv[DIM];
                            for (int m = 0; m < DIM; ++m)
                            {
Michael Turner's avatar
Michael Turner committed
356
                                basisDeriv[m] = *(m_derivUtils[typeIt->first]->VdmD[m])(k,typeIt->second[i]->NodeId(m_node->m_id));
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
                            }

                            for (int m = 0; m < DIM; ++m)
                            {
                                jacDetDeriv[m] = 0.0;
                                for (int n = 0; n < DIM; ++n)
                                {
                                    jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
                                }
                                jacDetDeriv[m] *= derivDet / fabs(typeIt->second[i]->maps[k][9]);
                            }

                            NekDouble jacDeriv[DIM][DIM][DIM];
                            for (int m = 0; m < DIM; ++m)
                            {
                                for (int n = 0; n < DIM; ++n)
                                {
                                    NekDouble delta = m == n ? 1.0 : 0.0;
                                    for (int l = 0; l < DIM; ++l)
                                    {
                                        jacDeriv[m][n][l] = delta * basisDeriv[l];
                                    }
                                }
                            }

382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
                            NekDouble jacDerivPhi[DIM][DIM][DIM];
                            for (int p = 0; p < DIM; ++p)
                            {
                                for (int m = 0; m < DIM; ++m)
                                {
                                    for (int n = 0; n < DIM; ++n)
                                    {
                                        jacDerivPhi[p][m][n] = 0.0;
                                        for (int l = 0; l < DIM; ++l)
                                        {
                                            // want phi_I^{-1} (l,n)
                                            jacDerivPhi[p][m][n] +=
                                                jacDeriv[p][m][l] * typeIt->second[i]->maps[k][l + 3*n];
                                        }
                                    }
                                }
                            }
399

400
401
                            NekDouble M2[DIM][DIM][DIM];
                            LEM2<DIM>(jacIdeal,jacDerivPhi,M2);
Michael Turner's avatar
Michael Turner committed
402

403
404
                            NekDouble M3[DIM][DIM][DIM][DIM];
                            LEM3<DIM>(jacDerivPhi,M3);
Michael Turner's avatar
Michael Turner committed
405
406

                            NekDouble frobProdA[DIM];
407
408
                            NekDouble frobProdB[DIM][DIM];
                            NekDouble frobProdC[DIM][DIM];
409
410
                            for (int m = 0; m < DIM; ++m)
                            {
411
412
413
414
415
416
                                frobProdA[m] = FrobProd<DIM>(M2[m],Emat);
                                for(int l = 0; l < DIM; l++)
                                {
                                    frobProdB[m][l] = FrobProd<DIM>(M3[m][l],Emat);
                                    frobProdC[m][l] = FrobProd<DIM>(M2[m],M2[l]);
                                }
417
418
419
420
                            }

                            for (int j = 0; j < DIM; ++j)
                            {
Michael Turner's avatar
Michael Turner committed
421
                                m_grad[j] += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) * (
Michael Turner's avatar
Michael Turner committed
422
                                        2.0 * mu * frobProdA[j] + K * lsigma * jacDetDeriv[j] / (2.0*sigma - jacDet));
423
424
                            }

425
                            int ct = 0;
Michael Turner's avatar
Michael Turner committed
426
                            for (int m = 0; m < DIM; ++m)
427
                            {
Michael Turner's avatar
Michael Turner committed
428
                                for(int l = m; l < DIM; ++l, ct++)
429
                                {
Michael Turner's avatar
Michael Turner committed
430
                                    m_grad[ct+DIM] += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) * (
431
                                        2.0 * mu * frobProdB[m][l] + 2.0 * mu * frobProdC[m][l] +
Michael Turner's avatar
Michael Turner committed
432
433
                                        jacDetDeriv[m]*jacDetDeriv[l]*K/(2.0*sigma-jacDet)/(2.0*sigma-jacDet)*(
                                            1.0 - jacDet*lsigma/(2.0*sigma-jacDet)));
434
                                }
435
                            }
436
                        }
437
                    }
Michael Turner's avatar
Michael Turner committed
438
439
440
441
442
443
444
                }
            }
            break;
        }

        case eHypEl:
        {
445
            const NekDouble nu = 0.45;
Michael Turner's avatar
Michael Turner committed
446
447
448
            const NekDouble mu = 1.0 / 2.0 / (1.0+nu);
            const NekDouble K  = 1.0 / 3.0 / (1.0 - 2.0 * nu);

Michael Turner's avatar
Michael Turner committed
449
            for(typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
Michael Turner's avatar
Michael Turner committed
450
            {
Michael Turner's avatar
Michael Turner committed
451
                NekVector<NekDouble> &quadW = m_derivUtils[typeIt->first]->quadW;
Michael Turner's avatar
Michael Turner committed
452
                for (int i = 0; i < typeIt->second.size(); ++i)
Michael Turner's avatar
Michael Turner committed
453
                {
Michael Turner's avatar
Michael Turner committed
454
                    for(int k = 0; k < m_derivUtils[typeIt->first]->pts; ++k)
455
                    {
Michael Turner's avatar
Michael Turner committed
456
                        jacDet = CalcIdealJac(i, k, derivs[typeIt->first], typeIt->second, jacIdeal);
457
                        minJacNew = min(minJacNew,jacDet);
Michael Turner's avatar
Michael Turner committed
458
                        NekDouble I1 = FrobeniusNorm(jacIdeal);
Michael Turner's avatar
Michael Turner committed
459

Michael Turner's avatar
Michael Turner committed
460
461
                        NekDouble sigma =
                            0.5*(jacDet + sqrt(jacDet*jacDet + 4.0*ep*ep));
462

463
464
465
466
                        if(sigma < numeric_limits<double>::min() && !gradient)
                        {
                            return numeric_limits<double>::max();
                        }
467
468
469
470
471
472

                        ASSERTL0(sigma > numeric_limits<double>::min(),
                                 std::string("dividing by zero ") +
                                 boost::lexical_cast<string>(sigma) + " " +
                                 boost::lexical_cast<string>(jacDet) + " " +
                                 boost::lexical_cast<string>(ep));
473

Michael Turner's avatar
Michael Turner committed
474
                        NekDouble lsigma = log(sigma);
Michael Turner's avatar
Michael Turner committed
475
                        integral += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) *
Michael Turner's avatar
Michael Turner committed
476
477
                                    (0.5 * mu * (I1 - 3.0 - 2.0*lsigma) +
                                     0.5 * K * lsigma * lsigma);
Michael Turner's avatar
Michael Turner committed
478

Michael Turner's avatar
Michael Turner committed
479
480
                        // Derivative of basis function in each direction
                        if(gradient)
Michael Turner's avatar
Michael Turner committed
481
                        {
Michael Turner's avatar
Michael Turner committed
482
483
484
485
486
                            NekDouble jacInvTrans[DIM][DIM];
                            NekDouble jacDetDeriv[DIM];

                            NekDouble phiM[DIM][DIM];
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
487
                            {
Michael Turner's avatar
Michael Turner committed
488
489
                                for (int n = 0; n < DIM; ++n)
                                {
Michael Turner's avatar
Michael Turner committed
490
                                    phiM[n][m] = derivs[typeIt->first][m][i][n][k];
Michael Turner's avatar
Michael Turner committed
491
                                }
Michael Turner's avatar
Michael Turner committed
492
493
                            }

Michael Turner's avatar
Michael Turner committed
494
495
                            InvTrans<DIM>(phiM, jacInvTrans);
                            NekDouble derivDet = Determinant<DIM>(phiM);
Michael Turner's avatar
Michael Turner committed
496

Michael Turner's avatar
Michael Turner committed
497
498
                            NekDouble basisDeriv[DIM];
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
499
                            {
Michael Turner's avatar
Michael Turner committed
500
                                basisDeriv[m] = *(m_derivUtils[typeIt->first]->VdmD[m])(k,typeIt->second[i]->NodeId(m_node->m_id));
Michael Turner's avatar
Michael Turner committed
501
502
                            }

Michael Turner's avatar
Michael Turner committed
503
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
504
                            {
Michael Turner's avatar
Michael Turner committed
505
506
                                jacDetDeriv[m] = 0.0;
                                for (int n = 0; n < DIM; ++n)
Michael Turner's avatar
Michael Turner committed
507
                                {
Michael Turner's avatar
Michael Turner committed
508
                                    jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
Michael Turner's avatar
Michael Turner committed
509
                                }
Michael Turner's avatar
Michael Turner committed
510
                                jacDetDeriv[m] *= derivDet / fabs(typeIt->second[i]->maps[k][9]);
Michael Turner's avatar
Michael Turner committed
511
512
                            }

Michael Turner's avatar
Michael Turner committed
513
                            NekDouble jacDeriv[DIM][DIM][DIM];
Michael Turner's avatar
Michael Turner committed
514
515
516
517
                            for (int m = 0; m < DIM; ++m)
                            {
                                for (int n = 0; n < DIM; ++n)
                                {
Michael Turner's avatar
Michael Turner committed
518
                                    NekDouble delta = m == n ? 1.0 : 0.0;
Michael Turner's avatar
Michael Turner committed
519
520
                                    for (int l = 0; l < DIM; ++l)
                                    {
Michael Turner's avatar
Michael Turner committed
521
                                        jacDeriv[m][n][l] = delta * basisDeriv[l];
Michael Turner's avatar
Michael Turner committed
522
523
524
525
                                    }
                                }
                            }

Michael Turner's avatar
Michael Turner committed
526
527
                            NekDouble jacDerivPhi[DIM][DIM][DIM];
                            for (int p = 0; p < DIM; ++p)
Michael Turner's avatar
Michael Turner committed
528
                            {
Michael Turner's avatar
Michael Turner committed
529
                                for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
530
                                {
Michael Turner's avatar
Michael Turner committed
531
532
533
534
535
536
537
538
539
540
                                    for (int n = 0; n < DIM; ++n)
                                    {
                                        jacDerivPhi[p][m][n] = 0.0;
                                        for (int l = 0; l < DIM; ++l)
                                        {
                                            // want phi_I^{-1} (l,n)
                                            jacDerivPhi[p][m][n] +=
                                                jacDeriv[p][m][l] * typeIt->second[i]->maps[k][l + 3*n];
                                        }
                                    }
Michael Turner's avatar
Michael Turner committed
541
542
543
                                }
                            }

Michael Turner's avatar
Michael Turner committed
544
                            NekDouble frobProd[DIM];
Michael Turner's avatar
Michael Turner committed
545
546
                            for (int m = 0; m < DIM; ++m)
                            {
Michael Turner's avatar
Michael Turner committed
547
548
549
550
551
                                frobProd[m] = FrobProd<DIM>(jacIdeal,jacDerivPhi[m]);
                            }

                            for (int j = 0; j < DIM; ++j)
                            {
Michael Turner's avatar
Michael Turner committed
552
                                m_grad[j] += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) * (
Michael Turner's avatar
Michael Turner committed
553
554
555
556
                                    mu * frobProd[j] + (jacDetDeriv[j] / (2.0*sigma - jacDet)
                                                        * (K * lsigma - mu)));
                            }

557
558
559

                            NekDouble frobProdHes[DIM][DIM]; //holder for the hessian frobprods
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
560
                            {
561
                                for(int l = m; l < DIM; ++l)
Michael Turner's avatar
Michael Turner committed
562
                                {
563
                                    frobProdHes[m][l] = FrobProd<DIM>(jacDerivPhi[m],jacDerivPhi[l]);
Michael Turner's avatar
Michael Turner committed
564
                                }
565
                            }
Michael Turner's avatar
Michael Turner committed
566

567
568
569
570
                            int ct = 0;
                            for (int m = 0; m < DIM; ++m)
                            {
                                for(int l = m; l < DIM; ++l, ct++)
Michael Turner's avatar
Michael Turner committed
571
                                {
Michael Turner's avatar
Michael Turner committed
572
                                    m_grad[ct+DIM] += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) * (
573
574
575
                                        mu * frobProdHes[m][l] +
                                        jacDetDeriv[m]*jacDetDeriv[l]/(2.0*sigma-jacDet)/(2.0*sigma-jacDet)*(
                                            K- jacDet*(K*lsigma-mu)/(2.0*sigma-jacDet)));
Michael Turner's avatar
Michael Turner committed
576
577
578
579
580
581
582
583
584
                                }
                            }
                        }
                    }
                }
            }
            break;
        }

Michael Turner's avatar
Michael Turner committed
585
        case eRoca:
Michael Turner's avatar
Michael Turner committed
586
        {
Michael Turner's avatar
Michael Turner committed
587
            for(typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
Michael Turner's avatar
Michael Turner committed
588
            {
Michael Turner's avatar
Michael Turner committed
589
                NekVector<NekDouble> &quadW = m_derivUtils[typeIt->first]->quadW;
Michael Turner's avatar
Michael Turner committed
590
                for (int i = 0; i < typeIt->second.size(); ++i)
Michael Turner's avatar
Michael Turner committed
591
                {
Michael Turner's avatar
Michael Turner committed
592
                    for(int k = 0; k < m_derivUtils[typeIt->first]->pts; ++k)
Michael Turner's avatar
Michael Turner committed
593
                    {
Michael Turner's avatar
Michael Turner committed
594
                        jacDet = CalcIdealJac(i, k, derivs[typeIt->first], typeIt->second, jacIdeal);
595
                        minJacNew = min(minJacNew,jacDet);
Michael Turner's avatar
Michael Turner committed
596
597
598
                        NekDouble frob = FrobeniusNorm(jacIdeal);
                        NekDouble sigma = 0.5*(jacDet +
                                        sqrt(jacDet*jacDet + 4.0*ep*ep));
Michael Turner's avatar
Michael Turner committed
599

600
601
602
603
                        if(sigma < numeric_limits<double>::min() && !gradient)
                        {
                            return numeric_limits<double>::max();
                        }
Michael Turner's avatar
Michael Turner committed
604
605
606
607
608
609
610

                        ASSERTL0(sigma > numeric_limits<double>::min(),
                                 std::string("dividing by zero ") +
                                 boost::lexical_cast<string>(sigma) + " " +
                                 boost::lexical_cast<string>(jacDet) + " " +
                                 boost::lexical_cast<string>(ep));

Michael Turner's avatar
Michael Turner committed
611
                        NekDouble W = frob / DIM / pow(fabs(sigma), 2.0/DIM);
Michael Turner's avatar
Michael Turner committed
612
                        integral += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) * W;
Michael Turner's avatar
Michael Turner committed
613

Michael Turner's avatar
Michael Turner committed
614
615
                        // Derivative of basis function in each direction
                        if(gradient)
Michael Turner's avatar
Michael Turner committed
616
                        {
Michael Turner's avatar
Michael Turner committed
617
618
619
620
621
                            NekDouble jacInvTrans[DIM][DIM];
                            NekDouble jacDetDeriv[DIM];

                            NekDouble phiM[DIM][DIM];
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
622
                            {
Michael Turner's avatar
Michael Turner committed
623
624
625
626
                                for (int n = 0; n < DIM; ++n)
                                {
                                    phiM[n][m] = derivs[typeIt->first][m][i][n][k];
                                }
Michael Turner's avatar
Michael Turner committed
627
628
                            }

Michael Turner's avatar
Michael Turner committed
629
630
                            InvTrans<DIM>(phiM, jacInvTrans);
                            NekDouble derivDet = Determinant<DIM>(phiM);
Michael Turner's avatar
Michael Turner committed
631

Michael Turner's avatar
Michael Turner committed
632
633
                            NekDouble basisDeriv[DIM];
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
634
                            {
Michael Turner's avatar
Michael Turner committed
635
                                basisDeriv[m] = *(m_derivUtils[typeIt->first]->VdmD[m])(k,typeIt->second[i]->NodeId(m_node->m_id));
Michael Turner's avatar
Michael Turner committed
636
637
                            }

Michael Turner's avatar
Michael Turner committed
638
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
639
                            {
Michael Turner's avatar
Michael Turner committed
640
641
                                jacDetDeriv[m] = 0.0;
                                for (int n = 0; n < DIM; ++n)
Michael Turner's avatar
Michael Turner committed
642
                                {
Michael Turner's avatar
Michael Turner committed
643
                                    jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
Michael Turner's avatar
Michael Turner committed
644
                                }
Michael Turner's avatar
Michael Turner committed
645
                                jacDetDeriv[m] *= derivDet / fabs(typeIt->second[i]->maps[k][9]);
Michael Turner's avatar
Michael Turner committed
646
647
                            }

Michael Turner's avatar
Michael Turner committed
648
                            NekDouble jacDeriv[DIM][DIM][DIM];
Michael Turner's avatar
Michael Turner committed
649
650
651
652
                            for (int m = 0; m < DIM; ++m)
                            {
                                for (int n = 0; n < DIM; ++n)
                                {
Michael Turner's avatar
Michael Turner committed
653
                                    NekDouble delta = m == n ? 1.0 : 0.0;
Michael Turner's avatar
Michael Turner committed
654
655
                                    for (int l = 0; l < DIM; ++l)
                                    {
Michael Turner's avatar
Michael Turner committed
656
                                        jacDeriv[m][n][l] = delta * basisDeriv[l];
Michael Turner's avatar
Michael Turner committed
657
658
659
660
                                    }
                                }
                            }

Michael Turner's avatar
Michael Turner committed
661
662
                            NekDouble jacDerivPhi[DIM][DIM][DIM];
                            for (int p = 0; p < DIM; ++p)
Michael Turner's avatar
Michael Turner committed
663
                            {
Michael Turner's avatar
Michael Turner committed
664
                                for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
665
                                {
Michael Turner's avatar
Michael Turner committed
666
667
668
669
670
671
672
673
674
675
                                    for (int n = 0; n < DIM; ++n)
                                    {
                                        jacDerivPhi[p][m][n] = 0.0;
                                        for (int l = 0; l < DIM; ++l)
                                        {
                                            // want phi_I^{-1} (l,n)
                                            jacDerivPhi[p][m][n] +=
                                                jacDeriv[p][m][l] * typeIt->second[i]->maps[k][l + 3*n];
                                        }
                                    }
Michael Turner's avatar
Michael Turner committed
676
677
678
                                }
                            }

Michael Turner's avatar
Michael Turner committed
679
                            NekDouble frobProd[DIM];
Michael Turner's avatar
Michael Turner committed
680
681
                            for (int m = 0; m < DIM; ++m)
                            {
Michael Turner's avatar
Michael Turner committed
682
683
684
685
686
                                frobProd[m] = FrobProd<DIM>(jacIdeal,jacDerivPhi[m]);
                            }

                            for (int j = 0; j < DIM; ++j)
                            {
Michael Turner's avatar
Michael Turner committed
687
688
                                m_grad[j] += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) * (
                                                2.0*W*(frobProd[j]/frob -
Michael Turner's avatar
Michael Turner committed
689
690
691
                                                jacDetDeriv[j]/DIM/(2.0*sigma-jacDet)));
                            }

692
693
                            NekDouble frobProdHes[DIM][DIM]; //holder for the hessian frobprods
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
694
                            {
695
                                for(int l = m; l < DIM; ++l)
Michael Turner's avatar
Michael Turner committed
696
                                {
697
                                    frobProdHes[m][l] = FrobProd<DIM>(jacDerivPhi[m],jacDerivPhi[l]);
Michael Turner's avatar
Michael Turner committed
698
                                }
699
                            }
Michael Turner's avatar
Michael Turner committed
700
701


702
703
704
705
                            int ct = 0;
                            for (int m = 0; m < DIM; ++m)
                            {
                                for(int l = m; l < DIM; ++l, ct++)
Michael Turner's avatar
Michael Turner committed
706
                                {
Michael Turner's avatar
Michael Turner committed
707
708
                                    m_grad[ct+DIM] += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) * (
                                            m_grad[m]*m_grad[l]/W + 2.0*W*(frobProdHes[m][l]/frob
709
710
                                            - 2.0 * frobProd[m]*frobProd[l]/frob/frob
                                            + jacDetDeriv[m]*jacDetDeriv[l] * jacDet/(2.0*sigma-jacDet)/(2.0*sigma-jacDet)/(2.0*sigma-jacDet)/DIM));
Michael Turner's avatar
Michael Turner committed
711
712
713
714
715
716
717
718
719
720
721
                                }
                            }
                        }
                    }
                }
            }
            break;
        }

        case eWins:
        {
Michael Turner's avatar
Michael Turner committed
722
            for(typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
Michael Turner's avatar
Michael Turner committed
723
            {
Michael Turner's avatar
Michael Turner committed
724
                NekVector<NekDouble> &quadW = m_derivUtils[typeIt->first]->quadW;
Michael Turner's avatar
Michael Turner committed
725
                for (int i = 0; i < typeIt->second.size(); ++i)
Michael Turner's avatar
Michael Turner committed
726
                {
Michael Turner's avatar
Michael Turner committed
727
                    for(int k = 0; k < m_derivUtils[typeIt->first]->pts; ++k)
Michael Turner's avatar
Michael Turner committed
728
                    {
Michael Turner's avatar
Michael Turner committed
729
                        jacDet = CalcIdealJac(i, k, derivs[typeIt->first], typeIt->second, jacIdeal);
730
                        minJacNew = min(minJacNew,jacDet);
Michael Turner's avatar
Michael Turner committed
731
732
733
                        NekDouble frob = FrobeniusNorm(jacIdeal);
                        NekDouble sigma = 0.5*(jacDet +
                                        sqrt(jacDet*jacDet + 4.0*ep*ep));
Michael Turner's avatar
Michael Turner committed
734

735
736
737
738
                        if(sigma < numeric_limits<double>::min() && !gradient)
                        {
                            return numeric_limits<double>::max();
                        }
Michael Turner's avatar
Michael Turner committed
739
740
741
742
743
744
745

                        ASSERTL0(sigma > numeric_limits<double>::min(),
                                 std::string("dividing by zero ") +
                                 boost::lexical_cast<string>(sigma) + " " +
                                 boost::lexical_cast<string>(jacDet) + " " +
                                 boost::lexical_cast<string>(ep));

Michael Turner's avatar
Michael Turner committed
746
                        NekDouble W = frob / sigma;
Michael Turner's avatar
Michael Turner committed
747
                        integral += quadW[k]* fabs(typeIt->second[i]->maps[k][9]) * W;
Michael Turner's avatar
Michael Turner committed
748

Michael Turner's avatar
Michael Turner committed
749
750
                        // Derivative of basis function in each direction
                        if(gradient)
Michael Turner's avatar
Michael Turner committed
751
                        {
Michael Turner's avatar
Michael Turner committed
752
753
754
755
756
                            NekDouble jacInvTrans[DIM][DIM];
                            NekDouble jacDetDeriv[DIM];

                            NekDouble phiM[DIM][DIM];
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
757
                            {
Michael Turner's avatar
Michael Turner committed
758
759
760
761
                                for (int n = 0; n < DIM; ++n)
                                {
                                    phiM[n][m] = derivs[typeIt->first][m][i][n][k];
                                }
Michael Turner's avatar
Michael Turner committed
762
763
                            }

Michael Turner's avatar
Michael Turner committed
764
765
                            InvTrans<DIM>(phiM, jacInvTrans);
                            NekDouble derivDet = Determinant<DIM>(phiM);
Michael Turner's avatar
Michael Turner committed
766

Michael Turner's avatar
Michael Turner committed
767
768
                            NekDouble basisDeriv[DIM];
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
769
                            {
Michael Turner's avatar
Michael Turner committed
770
                                basisDeriv[m] = *(m_derivUtils[typeIt->first]->VdmD[m])(k,typeIt->second[i]->NodeId(m_node->m_id));
Michael Turner's avatar
Michael Turner committed
771
772
                            }

Michael Turner's avatar
Michael Turner committed
773
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
774
                            {
Michael Turner's avatar
Michael Turner committed
775
776
                                jacDetDeriv[m] = 0.0;
                                for (int n = 0; n < DIM; ++n)
Michael Turner's avatar
Michael Turner committed
777
                                {
Michael Turner's avatar
Michael Turner committed
778
                                    jacDetDeriv[m] += jacInvTrans[m][n] * basisDeriv[n];
Michael Turner's avatar
Michael Turner committed
779
                                }
Michael Turner's avatar
Michael Turner committed
780
                                jacDetDeriv[m] *= derivDet / fabs(typeIt->second[i]->maps[k][9]);
Michael Turner's avatar
Michael Turner committed
781
782
                            }

Michael Turner's avatar
Michael Turner committed
783
                            NekDouble jacDeriv[DIM][DIM][DIM];
Michael Turner's avatar
Michael Turner committed
784
785
786
787
                            for (int m = 0; m < DIM; ++m)
                            {
                                for (int n = 0; n < DIM; ++n)
                                {
Michael Turner's avatar
Michael Turner committed
788
                                    NekDouble delta = m == n ? 1.0 : 0.0;
Michael Turner's avatar
Michael Turner committed
789
790
                                    for (int l = 0; l < DIM; ++l)
                                    {
Michael Turner's avatar
Michael Turner committed
791
                                        jacDeriv[m][n][l] = delta * basisDeriv[l];
Michael Turner's avatar
Michael Turner committed
792
793
794
795
                                    }
                                }
                            }

Michael Turner's avatar
Michael Turner committed
796
797
                            NekDouble jacDerivPhi[DIM][DIM][DIM];
                            for (int p = 0; p < DIM; ++p)
Michael Turner's avatar
Michael Turner committed
798
                            {
Michael Turner's avatar
Michael Turner committed
799
                                for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
800
                                {
Michael Turner's avatar
Michael Turner committed
801
802
803
804
805
806
807
808
809
810
                                    for (int n = 0; n < DIM; ++n)
                                    {
                                        jacDerivPhi[p][m][n] = 0.0;
                                        for (int l = 0; l < DIM; ++l)
                                        {
                                            // want phi_I^{-1} (l,n)
                                            jacDerivPhi[p][m][n] +=
                                                jacDeriv[p][m][l] * typeIt->second[i]->maps[k][l + 3*n];
                                        }
                                    }
Michael Turner's avatar
Michael Turner committed
811
812
813
                                }
                            }

Michael Turner's avatar
Michael Turner committed
814
                            NekDouble frobProd[DIM];
Michael Turner's avatar
Michael Turner committed
815
816
                            for (int m = 0; m < DIM; ++m)
                            {
Michael Turner's avatar
Michael Turner committed
817
818
819
820
821
                                frobProd[m] = FrobProd<DIM>(jacIdeal,jacDerivPhi[m]);
                            }

                            for (int j = 0; j < DIM; ++j)
                            {
Michael Turner's avatar
Michael Turner committed
822
                                m_grad[j] += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) * (
Michael Turner's avatar
Michael Turner committed
823
824
825
826
                                        W*(2.0*frobProd[j]/frob -
                                                jacDetDeriv[j]/(2.0*sigma-jacDet)));
                            }

827
828
                            NekDouble frobProdHes[DIM][DIM]; //holder for the hessian frobprods
                            for (int m = 0; m < DIM; ++m)
Michael Turner's avatar
Michael Turner committed
829
                            {
830
                                for(int l = m; l < DIM; ++l)
Michael Turner's avatar
Michael Turner committed
831
                                {
832
                                    frobProdHes[m][l] = FrobProd<DIM>(jacDerivPhi[m],jacDerivPhi[l]);
Michael Turner's avatar
Michael Turner committed
833
                                }
834
                            }
Michael Turner's avatar
Michael Turner committed
835
836


837
838
839
840
                            int ct = 0;
                            for (int m = 0; m < DIM; ++m)
                            {
                                for(int l = m; l < DIM; ++l, ct++)
Michael Turner's avatar
Michael Turner committed
841
                                {
Michael Turner's avatar
Michael Turner committed
842
843
                                    m_grad[ct+DIM] += quadW[k] * fabs(typeIt->second[i]->maps[k][9]) * (
                                        m_grad[m]*m_grad[l]/W + 2.0*W*(frobProdHes[m][l]/frob
844
845
                                            - 2.0 * frobProd[m]*frobProd[l]/frob/frob
                                            + 0.5*jacDetDeriv[m]*jacDetDeriv[l] * jacDet/(2.0*sigma-jacDet)/(2.0*sigma-jacDet)/(2.0*sigma-jacDet)));
Michael Turner's avatar
Michael Turner committed
846
847
848
849
850
851
852
                                }
                            }
                        }
                    }
                }
            }
            break;
Michael Turner's avatar
Michael Turner committed
853
        }
Michael Turner's avatar
Michael Turner committed
854
855
    }

Dave Moxey's avatar
Dave Moxey committed
856
    //ASSERTL0(std::isfinite(integral),"inf in integral");
857

Michael Turner's avatar
Michael Turner committed
858
    return integral;
Michael Turner's avatar
Michael Turner committed
859
    //return sqrt(m_grad[0]*m_grad[0] + m_grad[1]*m_grad[1]);
Michael Turner's avatar
Michael Turner committed
860
861
862
863
864
865
}

}
}

#endif