DiffusionLDGNS.cpp 32.9 KB
Newer Older
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
///////////////////////////////////////////////////////////////////////////////
//
// File: DiffusionLDGNS.cpp
//
// 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).
//
// 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: LDGNS diffusion class.
//
///////////////////////////////////////////////////////////////////////////////

#include "DiffusionLDGNS.h"
#include <iostream>
#include <iomanip>

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
39 40 41
#include <boost/core/ignore_unused.hpp>
#include <boost/algorithm/string/predicate.hpp>

42 43 44 45
#include <LocalRegions/Expansion2D.h>

namespace Nektar
{
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
46 47
std::string DiffusionLDGNS::type = GetDiffusionFactory().
RegisterCreatorFunction("LDGNS", DiffusionLDGNS::create);
48

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
49 50 51
DiffusionLDGNS::DiffusionLDGNS()
{
}
52

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
53 54 55 56 57 58
void DiffusionLDGNS::v_InitObject(
    LibUtilities::SessionReaderSharedPtr        pSession,
    Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
    m_session = pSession;
    m_session->LoadParameter ("Twall", m_Twall, 300.15);
59

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
60 61 62
    // Setting up the normals
    std::size_t nDim = pFields[0]->GetCoordim(0);
    std::size_t nTracePts = pFields[0]->GetTrace()->GetTotPoints();
63

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
64 65 66 67
    m_spaceDim = nDim;
    if (pSession->DefinesSolverInfo("HOMOGENEOUS"))
    {
        m_spaceDim = 3;
68 69
    }

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
70 71 72 73 74
    m_diffDim = m_spaceDim - nDim;

    m_traceVel = Array<OneD, Array<OneD, NekDouble> >{m_spaceDim};
    m_traceNormals = Array<OneD, Array<OneD, NekDouble> >{m_spaceDim};
    for(std::size_t i = 0; i < m_spaceDim; ++i)
75
    {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
76 77
        m_traceVel[i] = Array<OneD, NekDouble> {nTracePts, 0.0};
        m_traceNormals[i] = Array<OneD, NekDouble> {nTracePts};
78
    }
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
79
    pFields[0]->GetTrace()->GetNormals(m_traceNormals);
80

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
81 82 83 84 85 86
    // Create equation of state object
    std::string eosType;
    m_session->LoadSolverInfo("EquationOfState",
                              eosType, "IdealGas");
    m_eos = GetEquationOfStateFactory()
                            .CreateInstance(eosType, m_session);
87 88


Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
89 90 91 92
    // Set up {h} reference on the trace for penalty term
    //
    // Note, this shold be replaced with something smarter when merging
    // LDG with IP
93

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
94 95 96 97 98 99 100 101
    // Get min h per element
    std::size_t nElements = pFields[0]->GetExpSize();
    Array<OneD, NekDouble> hEle{nElements, 1.0};
    for (std::size_t e = 0; e < nElements; e++)
    {
        NekDouble h{1.0e+10};
        std::size_t expDim = pFields[0]->GetShapeDimension();
        switch(expDim)
102
        {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
103 104 105 106 107 108 109 110 111 112 113
            case 3:
            {
                LocalRegions::Expansion3DSharedPtr exp3D;
                exp3D = pFields[0]->GetExp(e)->as<LocalRegions::Expansion3D>();
                for(std::size_t i = 0; i < exp3D->GetNedges(); ++i)
                {
                    h = std::min(h, exp3D->GetGeom3D()->GetEdge(i)->GetVertex(0)->
                        dist(*(exp3D->GetGeom3D()->GetEdge(i)->GetVertex(1))));
                }
            break;
            }
114

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
115
            case 2:
116
            {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
117 118 119 120 121 122 123 124
                LocalRegions::Expansion2DSharedPtr exp2D;
                exp2D = pFields[0]->GetExp(e)->as<LocalRegions::Expansion2D>();
                for(std::size_t i = 0; i < exp2D->GetNedges(); ++i)
                {
                    h = std::min(h, exp2D->GetGeom2D()->GetEdge(i)->GetVertex(0)->
                        dist(*(exp2D->GetGeom2D()->GetEdge(i)->GetVertex(1))));
                }
            break;
125
            }
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
126 127 128 129
            case 1:
            {
                LocalRegions::Expansion1DSharedPtr exp1D;
                exp1D = pFields[0]->GetExp(e)->as<LocalRegions::Expansion1D>();
130

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
131 132
                h = std::min(h, exp1D->GetGeom1D()->GetVertex(0)->
                    dist(*(exp1D->GetGeom1D()->GetVertex(1))));
133

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
134 135 136
            break;
            }
            default:
137
            {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
138
                ASSERTL0(false,"Dimension out of bound.")
139
            }
Zhenguo Yan's avatar
Zhenguo Yan committed
140 141 142 143
            case 1:
            {
                LocalRegions::Expansion1DSharedPtr exp1D;
                exp1D = pFields[0]->GetExp(e)->as<LocalRegions::Expansion1D>();
144

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
        // Store scaling
        hEle[e] = h;
    }
    // Expand h from elements to points
    std::size_t nPts = pFields[0]->GetTotPoints();
    Array<OneD, NekDouble> hElePts{nPts, 0.0};
    Array<OneD, NekDouble> tmp;
    for (std::size_t e = 0; e < pFields[0]->GetExpSize(); e++)
    {
        std::size_t nElmtPoints     = pFields[0]->GetExp(e)->GetTotPoints();
        std::size_t physOffset      = pFields[0]->GetPhys_Offset(e);
        Vmath::Fill(nElmtPoints, hEle[e],
            tmp = hElePts + physOffset, 1);
    }
    // Get Fwd and Bwd traces
    Array<OneD, NekDouble> Fwd{nTracePts, 0.0};
    Array<OneD, NekDouble> Bwd{nTracePts, 0.0};
    pFields[0]->GetFwdBwdTracePhys(hElePts, Fwd, Bwd);
    // Fix boundaries
    std::size_t cnt = 0;
    std::size_t nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
    // Loop on the boundary regions
    for (std::size_t i = 0; i < nBndRegions; ++i)
    {
        // Number of boundary regions related to region 'i'
        std::size_t nBndEdges = pFields[0]->
        GetBndCondExpansions()[i]->GetExpSize();

        if (pFields[0]->GetBndConditions()[i]->GetBoundaryConditionType()
            == SpatialDomains::ePeriodic)
175
        {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
176
            continue;
177 178
        }

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
179 180
        // Get value from interior
        for (std::size_t e = 0; e < nBndEdges; ++e)
181
        {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
182 183
            std::size_t nBndEdgePts = pFields[0]->
            GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
184

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
185 186 187
            std::size_t id2 = pFields[0]->GetTrace()->
            GetPhys_Offset(pFields[0]->GetTraceMap()->
                           GetBndCondTraceToGlobalTraceMap(cnt++));
188

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
189
            Vmath::Vcopy(nBndEdgePts, &Fwd[id2], 1, &Bwd[id2], 1);
190 191
        }
    }
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213
    // Get average of traces
    Array<OneD, NekDouble> traceH{nTracePts, 1.0};
    m_traceOneOverH = Array<OneD, NekDouble> {nTracePts, 1.0};
    Vmath::Svtsvtp(nTracePts, 0.5, Fwd, 1, 0.5, Bwd, 1, traceH, 1);
    // Multiply by coefficient = - C11 / h
    m_session->LoadParameter ("LDGNSc11", m_C11, 1.0);
    Vmath::Sdiv(nTracePts, -m_C11, traceH, 1, m_traceOneOverH, 1);
}

/**
 * @brief Calculate weak DG Diffusion in the LDG form for the
 * Navier-Stokes (NS) equations:
 *
 * \f$ \langle\psi, \hat{u}\cdot n\rangle
 *   - \langle\nabla\psi \cdot u\rangle
 *     \langle\phi, \hat{q}\cdot n\rangle -
 *     (\nabla \phi \cdot q) \rangle \f$
 *
 * The equations that need a diffusion operator are those related
 * with the velocities and with the energy.
 *
 */
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
214
 void DiffusionLDGNS::v_Diffuse(
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
215 216 217
    const std::size_t                                  nConvectiveFields,
    const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
    const Array<OneD, Array<OneD, NekDouble> >        &inarray,
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
218
            Array<OneD, Array<OneD, NekDouble> >        &outarray,
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
219 220 221
    const Array<OneD, Array<OneD, NekDouble> >        &pFwd,
    const Array<OneD, Array<OneD, NekDouble> >        &pBwd)
{
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
222 223 224
    int nCoeffs   = fields[0]->GetNcoeffs();
    Array<OneD, Array<OneD, NekDouble> > tmp2(nConvectiveFields);
    for (int i = 0; i < nConvectiveFields; ++i)
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
225
    {
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
226
        tmp2[i] = Array<OneD, NekDouble>(nCoeffs, 0.0);
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
227
    }
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
228 229
    v_Diffuse_coeff(nConvectiveFields,fields,inarray,tmp2,pFwd,pBwd);
    for (int i = 0; i < nConvectiveFields; ++i)
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
230
    {
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
231
        fields[i]->BwdTrans             (tmp2[i], outarray[i]);
Yu Pan's avatar
Yu Pan committed
232
    }
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
233
}
Yu Pan's avatar
Yu Pan committed
234

Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257
// TODO:: REPLACED
void DiffusionLDGNS::v_Diffuse_coeff(
    const std::size_t                                  nConvectiveFields,
    const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
    const Array<OneD, Array<OneD, NekDouble> >        &inarray,
            Array<OneD, Array<OneD, NekDouble> >        &outarray,
    const Array<OneD, Array<OneD, NekDouble> >        &pFwd,
    const Array<OneD, Array<OneD, NekDouble> >        &pBwd)
{
    int i, j;
    int nDim      = fields[0]->GetCoordim(0);
    int nScalars  = inarray.num_elements();
    int nPts      = fields[0]->GetTotPoints();
    int nCoeffs   = fields[0]->GetNcoeffs();
    int nTracePts = fields[0]->GetTrace()->GetTotPoints();

    Array<OneD, NekDouble>               tmp1(nCoeffs);
    Array<OneD, Array<OneD, NekDouble> > tmp2(nConvectiveFields);

    Array<OneD, Array<OneD, Array<OneD, NekDouble> > > 
                                            derivativesO1(m_spaceDim);

    for (j = 0; j < m_spaceDim; ++j)
Yu Pan's avatar
Yu Pan committed
258
    {
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
259 260 261 262
        derivativesO1[j]   = Array<OneD, Array<OneD, NekDouble> >(
                                                            nScalars);

        for (i = 0; i < nScalars; ++i)
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
263
        {
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
264
            derivativesO1[j][i]   = Array<OneD, NekDouble>(nPts, 0.0);
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
265
        }
Yu Pan's avatar
Yu Pan committed
266 267
    }

Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
268 269
    DiffuseCalculateDerivative(fields,inarray,derivativesO1,pFwd,pBwd);
    
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
270 271
    // Initialisation viscous tensor
    m_viscTensor = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
272 273 274 275
                                                            (m_spaceDim);
    Array<OneD, Array<OneD, NekDouble> > viscousFlux(nConvectiveFields);

    for (j = 0; j < m_spaceDim; ++j)
Yu Pan's avatar
Yu Pan committed
276
    {
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
277 278 279
        m_viscTensor[j] = Array<OneD, Array<OneD, NekDouble> >(
                                                            nScalars+1);
        for (i = 0; i < nScalars+1; ++i)
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
280
        {
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
281
            m_viscTensor[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
282
        }
Yu Pan's avatar
Yu Pan committed
283
    }
Zhenguo Yan's avatar
Zhenguo Yan committed
284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318
    // Get average of traces
    Array<OneD, NekDouble> traceH{nTracePts, 1.0};
    m_traceOneOverH = Array<OneD, NekDouble> {nTracePts, 1.0};
    Vmath::Svtsvtp(nTracePts, 0.5, Fwd, 1, 0.5, Bwd, 1, traceH, 1);
    // Multiply by coefficient = - C11 / h
    m_session->LoadParameter ("LDGNSc11", m_C11, 1.0);
    Vmath::Sdiv(nTracePts, -m_C11, traceH, 1, m_traceOneOverH, 1);
}

/**
 * @brief Calculate weak DG Diffusion in the LDG form for the
 * Navier-Stokes (NS) equations:
 *
 * \f$ \langle\psi, \hat{u}\cdot n\rangle
 *   - \langle\nabla\psi \cdot u\rangle
 *     \langle\phi, \hat{q}\cdot n\rangle -
 *     (\nabla \phi \cdot q) \rangle \f$
 *
 * The equations that need a diffusion operator are those related
 * with the velocities and with the energy.
 *
 */
void DiffusionLDGNS::v_Diffuse(
    const std::size_t                                  nConvectiveFields,
    const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
    const Array<OneD, Array<OneD, NekDouble> >        &inarray,
          Array<OneD, Array<OneD, NekDouble> >        &outarray,
    const Array<OneD, Array<OneD, NekDouble> >        &pFwd,
    const Array<OneD, Array<OneD, NekDouble> >        &pBwd)
{
    std::size_t nDim      = fields[0]->GetCoordim(0);
    std::size_t nPts      = fields[0]->GetTotPoints();
    std::size_t nCoeffs   = fields[0]->GetNcoeffs();
    std::size_t nScalars  = inarray.num_elements();
    std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
Yu Pan's avatar
Yu Pan committed
319

Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
320
    for (i = 0; i < nConvectiveFields; ++i)
321
    {
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
322
        viscousFlux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
323
    }
324

Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
325
    DiffuseVolumeFlux(fields,inarray,derivativesO1,m_viscTensor);
326

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
327 328
    // Compute u from q_{\eta} and q_{\xi}
    // Obtain numerical fluxes
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
329 330
    DiffuseTraceFlux(fields,inarray,derivativesO1,m_viscTensor,viscousFlux,pFwd,pBwd);
    // v_NumericalFluxO2(fields, inarray, m_viscTensor, viscousFlux);
331

Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
332
    for (i = 0; i < nConvectiveFields; ++i)
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
333
    {
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
334
        tmp2[i] = Array<OneD, NekDouble>(nCoeffs, 0.0);
335

Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
336
        for (j = 0; j < nDim; ++j)
337
        {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
338 339
            fields[i]->IProductWRTDerivBase(j, m_viscTensor[j][i], tmp1);
            Vmath::Vadd(nCoeffs, tmp1, 1, tmp2[i], 1, tmp2[i], 1);
340
        }
Zhenguo Yan's avatar
Zhenguo Yan committed
341
    }
342

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
343 344 345 346
        // Evaulate  <\phi, \hat{F}\cdot n> - outarray[i]
        Vmath::Neg                      (nCoeffs, tmp2[i], 1);
        fields[i]->AddTraceIntegral     (viscousFlux[i], tmp2[i]);
        fields[i]->SetPhysState         (false);
Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
347 348 349 350 351 352 353 354 355 356 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 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399
        fields[i]->MultiplyByElmtInvMass(tmp2[i], outarray[i]);
    }
}

void DiffusionLDGNS::v_DiffuseCalculateDerivative(
    const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
    const Array<OneD, Array<OneD, NekDouble> >        &inarray,
    Array<OneD,Array<OneD, Array<OneD, NekDouble> > > &inarrayderivative,
    const Array<OneD, Array<OneD, NekDouble> >        &pFwd,
    const Array<OneD, Array<OneD, NekDouble> >        &pBwd)
{
    int nDim      = fields[0]->GetCoordim(0);
    int nScalars  = inarray.num_elements();
    int nCoeffs   = fields[0]->GetNcoeffs();
    int nTracePts = fields[0]->GetTrace()->GetTotPoints();
    int nConvectiveFields = fields.num_elements();

    Array<OneD, NekDouble>               tmp1(nCoeffs);
    Array<OneD, Array<OneD, NekDouble> > tmp2(nConvectiveFields);
    Array<OneD, Array<OneD, Array<OneD, NekDouble> > > numericalFluxO1(m_spaceDim);
                                                
    for (int j = 0; j < m_spaceDim; ++j)
    {
        numericalFluxO1[j] = Array<OneD, Array<OneD, NekDouble> >(nScalars);

        for (int i = 0; i < nScalars; ++i)
        {
            numericalFluxO1[j][i] = Array<OneD, NekDouble>(nTracePts, 0.0);
        }
    }

    NumericalFluxO1(fields, inarray, numericalFluxO1, pFwd, pBwd);

    for (std::size_t j = 0; j < nDim; ++j)
    {
        for (std::size_t i = 0; i < nScalars; ++i)
        {
            fields[i]->IProductWRTDerivBase (j, inarray[i], tmp1);
            Vmath::Neg                      (nCoeffs, tmp1, 1);
            fields[i]->AddTraceIntegral     (numericalFluxO1[j][i],
                                             tmp1);
            fields[i]->SetPhysState         (false);
            fields[i]->MultiplyByElmtInvMass(tmp1, tmp1);
            fields[i]->BwdTrans             (tmp1, inarrayderivative[j][i]);
        }
    }
    // For 3D Homogeneous 1D only take derivatives in 3rd direction
    if (m_diffDim == 1)
    {
        for (int i = 0; i < nScalars; ++i)
        {
            inarrayderivative[2][i] = m_homoDerivs[i];
        }
Zhenguo Yan's avatar
Zhenguo Yan committed
400 401 402 403 404 405 406

        // Evaulate  <\phi, \hat{F}\cdot n> - outarray[i]
        Vmath::Neg                      (nCoeffs, tmp2[i], 1);
        fields[i]->AddTraceIntegral     (viscousFlux[i], tmp2[i]);
        fields[i]->SetPhysState         (false);
        fields[i]->MultiplyByElmtInvMass(tmp2[i], tmp2[i]);
        fields[i]->BwdTrans             (tmp2[i], outarray[i]);
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
407 408 409
    }
}

Zhenguo Yan's avatar
merged  
Zhenguo Yan committed
410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436
void DiffusionLDGNS::v_DiffuseVolumeFlux(
    const Array<OneD, MultiRegions::ExpListSharedPtr>   &fields,
    const Array<OneD, Array<OneD, NekDouble>>           &inarray,
    Array<OneD,Array<OneD, Array<OneD, NekDouble> > >   &inarrayderivative,
    Array<OneD, Array<OneD, Array<OneD, NekDouble> > >  &VolumeFlux,
    Array< OneD, int >                                  &nonZeroIndex) 
{

    boost::ignore_unused(fields,nonZeroIndex);
    m_fluxVectorNS(inarray, inarrayderivative, VolumeFlux);

}

void DiffusionLDGNS::v_DiffuseTraceFlux(
    const Array<OneD, MultiRegions::ExpListSharedPtr>   &fields,
    const Array<OneD, Array<OneD, NekDouble>>           &inarray,
    Array<OneD,Array<OneD, Array<OneD, NekDouble> > >   &inarrayderivative,
    Array<OneD, Array<OneD, Array<OneD, NekDouble> > >  &VolumeFlux,
    Array<OneD, Array<OneD, NekDouble> >                &TraceFlux,
    const Array<OneD, Array<OneD, NekDouble>>           &pFwd,
    const Array<OneD, Array<OneD, NekDouble>>           &pBwd,
    Array< OneD, int >                                  &nonZeroIndex)
{
    boost::ignore_unused(inarray,inarrayderivative,nonZeroIndex);
    NumericalFluxO2(fields, VolumeFlux, TraceFlux, pFwd, pBwd);
}

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
437 438 439 440 441 442 443 444 445 446 447 448 449 450
/**
 * @brief Builds the numerical flux for the 1st order derivatives
 *
 */
void DiffusionLDGNS::NumericalFluxO1(
    const Array<OneD, MultiRegions::ExpListSharedPtr>        &fields,
    const Array<OneD, Array<OneD, NekDouble> >               &inarray,
          Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
                                                    &numericalFluxO1,
    const Array<OneD, Array<OneD, NekDouble> >               &pFwd,
    const Array<OneD, Array<OneD, NekDouble> >               &pBwd)
{
    std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
    std::size_t nScalars  = inarray.num_elements();
451

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
452 453 454 455 456 457
    //Upwind
    Array<OneD, Array<OneD, NekDouble> > numflux{nScalars};
    for (std::size_t i = 0; i < nScalars; ++i)
    {
        numflux[i] = {pFwd[i]};
    }
458

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
459 460 461 462 463
    // Modify the values in case of boundary interfaces
    if (fields[0]->GetBndCondExpansions().num_elements())
    {
        ApplyBCsO1(fields, inarray, pFwd, pBwd, numflux);
    }
464

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
465 466 467 468
    // Splitting the numerical flux into the dimensions
    for (std::size_t j = 0; j < m_spaceDim; ++j)
    {
        for (std::size_t i = 0; i < nScalars; ++i)
469
        {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
470 471
            Vmath::Vmul(nTracePts, m_traceNormals[j], 1,numflux[i], 1,
                numericalFluxO1[j][i], 1);
472 473
        }
    }
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
474 475 476 477 478 479 480 481 482 483 484 485 486 487
}

/**
 * @brief Imposes appropriate bcs for the 1st order derivatives
 *
 */
void DiffusionLDGNS::ApplyBCsO1(
    const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
    const Array<OneD, Array<OneD, NekDouble> >        &inarray,
    const Array<OneD, Array<OneD, NekDouble> >        &pFwd,
    const Array<OneD, Array<OneD, NekDouble> >        &pBwd,
          Array<OneD, Array<OneD, NekDouble> >        &fluxO1)
{
    boost::ignore_unused(pBwd);
488

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
489 490
    std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
    std::size_t nScalars  = inarray.num_elements();
491

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510
    Array<OneD, NekDouble> tmp1{nTracePts, 0.0};
    Array<OneD, NekDouble> tmp2{nTracePts, 0.0};
    Array<OneD, NekDouble> Tw{nTracePts, m_Twall};

    Array< OneD, Array<OneD, NekDouble > > scalarVariables{nScalars};

    // Extract internal values of the scalar variables for Neumann bcs
    for (std::size_t i = 0; i < nScalars; ++i)
    {
        scalarVariables[i] = Array<OneD, NekDouble>{nTracePts, 0.0};
    }

    // Compute boundary conditions for velocities
    for (std::size_t i = 0; i < nScalars-1; ++i)
    {
        // Note that cnt has to loop on nBndRegions and nBndEdges
        // and has to be reset to zero per each equation
        std::size_t cnt = 0;
        std::size_t nBndRegions = fields[i+1]->
511
            GetBndCondExpansions().num_elements();
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
512 513 514 515
        for (std::size_t j = 0; j < nBndRegions; ++j)
        {
            if (fields[i+1]->GetBndConditions()[j]->
                GetBoundaryConditionType() == SpatialDomains::ePeriodic)
516
            {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
517 518
                continue;
            }
519

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
520
            std::size_t nBndEdges = fields[i+1]->
521
                GetBndCondExpansions()[j]->GetExpSize();
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
522 523 524
            for (std::size_t e = 0; e < nBndEdges; ++e)
            {
                std::size_t nBndEdgePts = fields[i+1]->
525 526
                    GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
527
                std::size_t id1 = fields[i+1]->
528 529
                    GetBndCondExpansions()[j]->GetPhys_Offset(e);

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
530
                std::size_t id2 = fields[0]->GetTrace()->
531
                    GetPhys_Offset(fields[0]->GetTraceMap()->
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
532
                    GetBndCondTraceToGlobalTraceMap(cnt++));
533

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
534 535 536 537 538 539 540
                if (boost::iequals(fields[i]->GetBndConditions()[j]->
                    GetUserDefined(),"WallViscous") ||
                    boost::iequals(fields[i]->GetBndConditions()[j]->
                    GetUserDefined(),"WallAdiabatic"))
                {
                    // Reinforcing bcs for velocity in case of Wall bcs
                    Vmath::Zero(nBndEdgePts, &scalarVariables[i][id2], 1);
541

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
542 543 544 545 546 547 548 549 550 551
                }
                else if (
                    boost::iequals(fields[i]->GetBndConditions()[j]->
                    GetUserDefined(),"Wall") ||
                    boost::iequals(fields[i]->GetBndConditions()[j]->
                    GetUserDefined(),"Symmetry"))
                {
                    // Symmetry bc: normal velocity is zero
                    //    get all velocities at once because we need u.n
                    if (i==0)
552
                    {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
553 554 555
                        //    tmp1 = -(u.n)
                        Vmath::Zero(nBndEdgePts, tmp1, 1);
                        for (std::size_t k = 0; k < nScalars-1; ++k)
556
                        {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
557 558
                            Vmath::Vdiv(nBndEdgePts,
                              &(fields[k+1]->GetBndCondExpansions()[j]->
559
                                  UpdatePhys())[id1], 1,
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
560
                              &(fields[0]->GetBndCondExpansions()[j]->
561
                                  UpdatePhys())[id1], 1,
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
562 563 564 565 566 567 568 569 570 571
                              &scalarVariables[k][id2], 1);
                            Vmath::Vvtvp(nBndEdgePts,
                                    &m_traceNormals[k][id2], 1,
                                    &scalarVariables[k][id2], 1,
                                    &tmp1[0], 1,
                                    &tmp1[0], 1);
                        }
                        Vmath::Smul(nBndEdgePts, -1.0,
                                    &tmp1[0], 1,
                                    &tmp1[0], 1);
572

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
573 574 575 576 577 578 579 580 581
                        //    u_i - (u.n)n_i
                        for (std::size_t k = 0; k < nScalars-1; ++k)
                        {
                            Vmath::Vvtvp(nBndEdgePts,
                                    &tmp1[0], 1,
                                    &m_traceNormals[k][id2], 1,
                                    &scalarVariables[k][id2], 1,
                                    &scalarVariables[k][id2], 1);
                        }
582 583 584
                    }
                }
                else if (fields[i]->GetBndConditions()[j]->
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
585
                         GetBoundaryConditionType() ==
586 587
                         SpatialDomains::eDirichlet)
                {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
588 589 590 591 592 593 594
                    // Imposing velocity bcs if not Wall
                    Vmath::Vdiv(nBndEdgePts,
                            &(fields[i+1]->GetBndCondExpansions()[j]->
                              UpdatePhys())[id1], 1,
                            &(fields[0]->GetBndCondExpansions()[j]->
                              UpdatePhys())[id1], 1,
                            &scalarVariables[i][id2], 1);
595 596 597
                }

                // For Dirichlet boundary condition: uflux = u_bcs
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
598
                if (fields[i]->GetBndConditions()[j]->
599
                    GetBoundaryConditionType() ==
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
600
                    SpatialDomains::eDirichlet)
601
                {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
602 603 604
                    Vmath::Vcopy(nBndEdgePts,
                                 &scalarVariables[i][id2], 1,
                                 &fluxO1[i][id2], 1);
605 606 607
                }

                // For Neumann boundary condition: uflux = u_+
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
608 609 610
                else if ((fields[i]->GetBndConditions()[j])->
                         GetBoundaryConditionType() ==
                         SpatialDomains::eNeumann)
611
                {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
612 613 614
                    Vmath::Vcopy(nBndEdgePts,
                                 &pFwd[i][id2], 1,
                                 &fluxO1[i][id2], 1);
615
                }
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
616 617 618 619 620 621 622 623 624 625 626 627 628 629 630

                // Building kinetic energy to be used for T bcs
                Vmath::Vmul(nBndEdgePts,
                            &scalarVariables[i][id2], 1,
                            &scalarVariables[i][id2], 1,
                            &tmp1[id2], 1);

                Vmath::Smul(nBndEdgePts, 0.5,
                            &tmp1[id2], 1,
                            &tmp1[id2], 1);

                Vmath::Vadd(nBndEdgePts,
                            &tmp2[id2], 1,
                            &tmp1[id2], 1,
                            &tmp2[id2], 1);
631 632 633 634
            }
        }
    }

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
635 636 637 638 639
    // Compute boundary conditions  for temperature
    std::size_t cnt = 0;
    std::size_t nBndRegions = fields[nScalars]->
    GetBndCondExpansions().num_elements();
    for (std::size_t j = 0; j < nBndRegions; ++j)
640
    {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
641 642 643
        if (fields[nScalars]->GetBndConditions()[j]->
            GetBoundaryConditionType() ==
            SpatialDomains::ePeriodic)
644
        {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
645
            continue;
646 647
        }

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
648 649 650
        std::size_t nBndEdges = fields[nScalars]->
        GetBndCondExpansions()[j]->GetExpSize();
        for (std::size_t e = 0; e < nBndEdges; ++e)
651
        {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
652 653
            std::size_t nBndEdgePts = fields[nScalars]->
            GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();
654

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
655 656
            std::size_t id1 = fields[nScalars]->
            GetBndCondExpansions()[j]->GetPhys_Offset(e);
657

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
658 659 660
            std::size_t id2 = fields[0]->GetTrace()->
            GetPhys_Offset(fields[0]->GetTraceMap()->
                           GetBndCondTraceToGlobalTraceMap(cnt++));
661

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678
            // Imposing Temperature Twall at the wall
            if (boost::iequals(fields[nScalars]->GetBndConditions()[j]->
                GetUserDefined(),"WallViscous"))
            {
                Vmath::Vcopy(nBndEdgePts,
                             &Tw[0], 1,
                             &scalarVariables[nScalars-1][id2], 1);
            }
            // Imposing Temperature through condition on the Energy
            // for no wall boundaries (e.g. farfield)
            else if (fields[nScalars]->GetBndConditions()[j]->
                     GetBoundaryConditionType() ==
                     SpatialDomains::eDirichlet)
            {
                // Use equation of state to evaluate temperature
                NekDouble rho, ene;
                for(std::size_t n = 0; n < nBndEdgePts; ++n)
679
                {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
680 681 682 683 684 685 686 687
                    rho = fields[0]->
                              GetBndCondExpansions()[j]->
                              GetPhys()[id1+n];
                    ene = fields[nScalars]->
                              GetBndCondExpansions()[j]->
                              GetPhys()[id1 +n] / rho - tmp2[id2+n];
                    scalarVariables[nScalars-1][id2+n] =
                            m_eos->GetTemperature(rho, ene);
688
                }
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711
            }

            // For Dirichlet boundary condition: uflux = u_bcs
            if (fields[nScalars]->GetBndConditions()[j]->
                GetBoundaryConditionType() == SpatialDomains::eDirichlet &&
                !boost::iequals(fields[nScalars]->GetBndConditions()[j]
                ->GetUserDefined(), "WallAdiabatic"))
            {
                Vmath::Vcopy(nBndEdgePts,
                             &scalarVariables[nScalars-1][id2], 1,
                             &fluxO1[nScalars-1][id2], 1);

            }

            // For Neumann boundary condition: uflux = u_+
            else if (((fields[nScalars]->GetBndConditions()[j])->
                      GetBoundaryConditionType() == SpatialDomains::eNeumann) ||
                      boost::iequals(fields[nScalars]->GetBndConditions()[j]->
                                    GetUserDefined(), "WallAdiabatic"))
            {
                Vmath::Vcopy(nBndEdgePts,
                             &pFwd[nScalars-1][id2], 1,
                             &fluxO1[nScalars-1][id2], 1);
712 713 714 715

            }
        }
    }
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
716 717 718 719 720 721 722 723 724 725 726 727 728 729 730
}

/**
 * @brief Build the numerical flux for the 2nd order derivatives
 *
 */
void DiffusionLDGNS::NumericalFluxO2(
    const Array<OneD, MultiRegions::ExpListSharedPtr>        &fields,
          Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &qfield,
          Array<OneD, Array<OneD, NekDouble> >               &qflux,
    const Array<OneD, Array<OneD, NekDouble> >               &uFwd,
    const Array<OneD, Array<OneD, NekDouble> >               &uBwd)
{
    std::size_t nTracePts  = fields[0]->GetTrace()->GetTotPoints();
    std::size_t nVariables = fields.num_elements();
731

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
732 733 734
    // Initialize penalty flux
    Array<OneD, Array<OneD, NekDouble> >  fluxPen{nVariables-1};
    for (std::size_t i = 0; i < nVariables-1; ++i)
735
    {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
736 737 738 739 740
        fluxPen[i] = Array<OneD, NekDouble>{nTracePts, 0.0};
    }

    // Get penalty flux
    m_fluxPenaltyNS(uFwd, uBwd, fluxPen);
741

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
742 743 744
    // Evaluate Riemann flux
    // qflux = \hat{q} \cdot u = q \cdot n
    // Notice: i = 1 (first row of the viscous tensor is zero)
745

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
746 747 748 749 750 751 752 753 754
    Array<OneD, NekDouble > qFwd{nTracePts};
    Array<OneD, NekDouble > qBwd{nTracePts};
    Array<OneD, NekDouble > qtemp{nTracePts};
    Array<OneD, NekDouble > qfluxtemp{nTracePts, 0.0};
    std::size_t nDim = fields[0]->GetCoordim(0);
    for (std::size_t i = 1; i < nVariables; ++i)
    {
        qflux[i] = Array<OneD, NekDouble> {nTracePts, 0.0};
        for (std::size_t j = 0; j < nDim; ++j)
755
        {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
756 757 758 759 760 761 762 763 764
            // Compute qFwd and qBwd value of qfield in position 'ji'
            fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);

            // Downwind
            Vmath::Vcopy(nTracePts, qBwd, 1, qfluxtemp, 1);

            // Multiply the Riemann flux by the trace normals
            Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
                        qfluxtemp, 1);
765

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
766 767 768 769 770 771
            // Add penalty term
            Vmath::Vvtvp(nTracePts, m_traceOneOverH, 1, fluxPen[i-1], 1,
                qfluxtemp, 1, qfluxtemp, 1);

            // Impose weak boundary condition with flux
            if (fields[0]->GetBndCondExpansions().num_elements())
772
            {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
773
                ApplyBCsO2(fields, i, j, qfield[j][i], qFwd, qBwd, qfluxtemp);
774 775
            }

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828
            // Store the final flux into qflux
            Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
        }
    }
}


/**
 * @brief Imposes appropriate bcs for the 2nd order derivatives
 *
 */
void DiffusionLDGNS::ApplyBCsO2(
    const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
    const std::size_t                                  var,
    const std::size_t                                  dir,
    const Array<OneD, const NekDouble>                &qfield,
    const Array<OneD, const NekDouble>                &qFwd,
    const Array<OneD, const NekDouble>                &qBwd,
          Array<OneD,       NekDouble>                &penaltyflux)
{
    boost::ignore_unused(qfield, qBwd);

    std::size_t cnt = 0;
    std::size_t nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
    // Loop on the boundary regions to apply appropriate bcs
    for (std::size_t i = 0; i < nBndRegions; ++i)
    {
        // Number of boundary regions related to region 'i'
        std::size_t nBndEdges = fields[var]->
        GetBndCondExpansions()[i]->GetExpSize();

        if (fields[var]->GetBndConditions()[i]->GetBoundaryConditionType()
            == SpatialDomains::ePeriodic)
        {
            continue;
        }

        // Weakly impose bcs by modifying flux values
        for (std::size_t e = 0; e < nBndEdges; ++e)
        {
            std::size_t nBndEdgePts = fields[var]->
            GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();

            std::size_t id2 = fields[0]->GetTrace()->
            GetPhys_Offset(fields[0]->GetTraceMap()->
                           GetBndCondTraceToGlobalTraceMap(cnt++));

            // In case of Dirichlet bcs:
            // uflux = gD
            if (fields[var]->GetBndConditions()[i]->
               GetBoundaryConditionType() == SpatialDomains::eDirichlet
               && !boost::iequals(fields[var]->GetBndConditions()[i]->
                                  GetUserDefined(), "WallAdiabatic"))
829
            {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855
                Vmath::Vmul(nBndEdgePts,
                            &m_traceNormals[dir][id2], 1,
                            &qFwd[id2], 1,
                            &penaltyflux[id2], 1);
            }
            // 3.4) In case of Neumann bcs:
            // uflux = u+
            else if ((fields[var]->GetBndConditions()[i])->
                GetBoundaryConditionType() == SpatialDomains::eNeumann)
            {
                ASSERTL0(false,
                         "Neumann bcs not implemented for LDGNS");

                /*
                Vmath::Vmul(nBndEdgePts,
                            &m_traceNormals[dir][id2], 1,
                            &(fields[var]->
                              GetBndCondExpansions()[i]->
                              UpdatePhys())[id1], 1,
                            &penaltyflux[id2], 1);
                 */
            }
            else if (boost::iequals(fields[var]->GetBndConditions()[i]->
                                   GetUserDefined(), "WallAdiabatic"))
            {
                if ((var == m_spaceDim + 1))
856
                {
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
857
                    Vmath::Zero(nBndEdgePts, &penaltyflux[id2], 1);
858
                }
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
859
                else
860 861
                {

Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
862 863 864 865
                    Vmath::Vmul(nBndEdgePts,
                                &m_traceNormals[dir][id2], 1,
                                &qFwd[id2], 1,
                                &penaltyflux[id2], 1);
866 867 868 869 870

                }
            }
        }
    }
Zhenguo Yan's avatar
TMP  
Zhenguo Yan committed
871
}
872
}//end of namespace Nektar