From 76ba87a6c2fafac42c64d4e2367dda252bc8addd Mon Sep 17 00:00:00 2001 From: Chris Cantwell Date: Mon, 10 Feb 2014 18:06:16 +0000 Subject: [PATCH] Cleaned up formatting. --- .../FieldConvert/ProcessQCriterion.cpp | 299 ++++++++++-------- 1 file changed, 165 insertions(+), 134 deletions(-) diff --git a/utilities/PostProcessing/FieldConvert/ProcessQCriterion.cpp b/utilities/PostProcessing/FieldConvert/ProcessQCriterion.cpp index 28be24fd9..ac682a18d 100644 --- a/utilities/PostProcessing/FieldConvert/ProcessQCriterion.cpp +++ b/utilities/PostProcessing/FieldConvert/ProcessQCriterion.cpp @@ -47,25 +47,26 @@ namespace Nektar namespace Utilities { ModuleKey ProcessQCriterion::className = - GetModuleFactory().RegisterCreatorFunction( - ModuleKey(eProcessModule, "QCriterion"), - ProcessQCriterion::create, "Computes Q-Criterion."); - - ProcessQCriterion::ProcessQCriterion(FieldSharedPtr f) : ProcessModule(f) + GetModuleFactory().RegisterCreatorFunction( + ModuleKey(eProcessModule, "QCriterion"), + ProcessQCriterion::create, "Computes Q-Criterion."); + + ProcessQCriterion::ProcessQCriterion(FieldSharedPtr f) : + ProcessModule(f) { } - + ProcessQCriterion::~ProcessQCriterion() { } - + void ProcessQCriterion::Process(po::variables_map &vm) { if (m_f->m_verbose) { cout << "ProcessQCriterion: Calculating Q Criterion..." << endl; } - + int i, j; int expdim = m_f->m_graph->GetMeshDimension(); int spacedim = expdim; @@ -77,146 +78,177 @@ namespace Nektar int nfields = m_f->m_fielddef[0]->m_fields.size(); if (spacedim == 1 || spacedim == 2) { - cerr << "\n Error: ProcessQCriterion must be computed for a 3D (or quasi-3D) case. \n" << endl; + cerr << "\n Error: ProcessQCriterion must be computed for a 3D" + " (or quasi-3D) case. \n" << endl; } - int addfields = 1; //For calculating Q-Criterion only 1 field must be added - + + //For calculating Q-Criterion only 1 field must be added + int addfields = 1; + int npoints = m_f->m_exp[0]->GetNpoints(); - - Array > grad(nfields * nfields); - - Array > omega(nfields * nfields); - Array > S (nfields * nfields); - - Array > outfield (addfields); - Array > outfield1(addfields); - Array > outfield2(addfields); - Array > outfield3(addfields); - - + + Array > grad(nfields * nfields); + + Array > omega(nfields * nfields); + Array > S (nfields * nfields); + + Array > outfield (addfields); + Array > outfield1(addfields); + Array > outfield2(addfields); + Array > outfield3(addfields); + m_f->m_exp.resize(nfields+addfields); - - + for (i = 0; i < nfields*nfields; ++i) { grad[i] = Array(npoints); } - + for (i = 0; i < addfields; ++i) { - outfield[i] = Array(npoints); //Will store the Q-Criterion - outfield1[i] = Array(npoints); - outfield2[i] = Array(npoints); - outfield3[i] = Array(npoints); - - omega[i] = Array(npoints); - S[i] = Array(npoints); + //Will store the Q-Criterion + outfield[i] = Array(npoints); + outfield1[i] = Array(npoints); + outfield2[i] = Array(npoints); + outfield3[i] = Array(npoints); + + omega[i] = Array(npoints); + S[i] = Array(npoints); + } + + for (i = 0; i < nfields; ++i) + { + m_f->m_exp[i]->PhysDeriv(m_f->m_exp[i]->GetPhys(), + grad[i*nfields], + grad[i*nfields+1], + grad[i*nfields+2]); + } + + // W_x = Wy - Vz + Vmath::Vsub(npoints, grad[2 * nfields + 1], 1, + grad[1 * nfields + 2], 1, + outfield1[0], 1); + // W_x^2 + Vmath::Vmul(npoints, outfield1[0], 1, + outfield1[0], 1, + outfield1[0], 1); + + // W_y = Uz - Wx + Vmath::Vsub(npoints, grad[0 * nfields + 2], 1, + grad[2 * nfields + 0], 1, + outfield2[0], 1); + // W_y^2 + Vmath::Vmul(npoints, outfield2[0], 1, + outfield2[0], 1, + outfield2[0], 1); + + // W_z = Vx - Uy + Vmath::Vsub(npoints, grad[1 * nfields + 0], 1, + grad[0 * nfields + 1], 1, + outfield3[0], 1); + // W_z^2 + Vmath::Vmul(npoints, outfield3[0], 1, + outfield3[0], 1, + outfield3[0], 1); + + // add fields omega = 0.5*(W_x^2 + W_y^2 + W_z^2) + + NekDouble fac = 0.5; + Vmath::Vadd(npoints, &outfield1[0][0], 1, + &outfield2[0][0], 1, + &omega[0][0], 1); + Vmath::Vadd(npoints, &omega[0][0], 1, + &outfield3[0][0], 1, + &omega[0][0], 1); + + for (int k = 0; k < addfields; ++k) + { + Vmath::Smul(npoints, fac, &omega[k][0], 1, &omega[k][0], 1); + } + + Vmath::Zero(npoints, &outfield1[0][0], 1); + Vmath::Zero(npoints, &outfield2[0][0], 1); + Vmath::Zero(npoints, &outfield3[0][0], 1); + + Vmath::Vmul(npoints, grad[0 * nfields + 0], 1, + grad[0 * nfields + 0], 1, + outfield1[0], 1); + Vmath::Vmul(npoints, grad[1 * nfields + 1], 1, + grad[1 * nfields + 1], 1, + outfield2[0], 1); + Vmath::Vmul(npoints, grad[2 * nfields + 2], 1, + grad[2 * nfields + 2], 1, + outfield3[0], 1); + + Vmath::Vadd(npoints, &outfield1[0][0], 1, + &outfield2[0][0], 1, + &S[0][0], 1); + Vmath::Vadd(npoints, &S[0][0], 1, + &outfield3[0][0], 1, + &S[0][0], 1); + + // W_y + V_z + Vmath::Vadd(npoints, grad[2 * nfields + 1], 1, + grad[1 * nfields + 2], 1, + outfield1[0], 1); + Vmath::Vmul(npoints, &outfield1[0][0], 1, + &outfield1[0][0], 1, + &outfield1[0][0], 1); + + // U_z + W_x + Vmath::Vadd(npoints, grad[0 * nfields + 2], 1, + grad[2 * nfields + 0], 1, + outfield2[0], 1); + Vmath::Vmul(npoints, &outfield2[0][0], 1, + &outfield2[0][0], 1, + &outfield2[0][0], 1); + + // V_x + U_y + Vmath::Vadd(npoints, grad[1 * nfields + 0], 1, + grad[0 * nfields + 1], 1, + outfield3[0], 1); + Vmath::Vmul(npoints, &outfield3[0][0], 1, + &outfield3[0][0], 1, + &outfield3[0][0], 1); + + Vmath::Vadd(npoints, &outfield1[0][0], 1, + &outfield2[0][0], 1, + &outfield2[0][0], 1); + Vmath::Vadd(npoints, &outfield2[0][0], 1, + &outfield3[0][0], 1, + &outfield3[0][0], 1); + + for (int k = 0; k < addfields; ++k) + { + Vmath::Smul(npoints, fac, &outfield3[k][0], 1, + &outfield3[k][0], 1); + } + + Vmath::Vadd(npoints, &outfield3[0][0], 1, &S[0][0], 1, &S[0][0], 1); + Vmath::Vsub(npoints, omega[0], 1, S[0], 1, outfield[0], 1); + + for (int k = 0; k < addfields; ++k) + { + Vmath::Smul(npoints, fac, &outfield[k][0], 1, + &outfield[k][0], 1); } - - for (i = 0; i < nfields; ++i) - { - - m_f->m_exp[i]->PhysDeriv(m_f->m_exp[i]->GetPhys(), - grad[i*nfields], - grad[i*nfields+1], - grad[i*nfields+2]); - } - - // W_x = Wy - Vz - Vmath::Vsub(npoints, grad[2 * nfields + 1], 1, grad[1 * nfields + 2], 1, - outfield1[0], 1); - // W_x^2 - Vmath::Vmul(npoints, outfield1[0], 1, outfield1[0], 1, outfield1[0], 1); - - // W_y = Uz - Wx - Vmath::Vsub(npoints, grad[0 * nfields + 2], 1, grad[2 * nfields + 0], 1, - outfield2[0], 1); - // W_y^2 - Vmath::Vmul(npoints, outfield2[0], 1, outfield2[0], 1, outfield2[0], 1); - - // W_z = Vx - Uy - Vmath::Vsub(npoints, grad[1 * nfields + 0], 1, grad[0 * nfields + 1], 1, - outfield3[0], 1); - // W_z^2 - Vmath::Vmul(npoints, outfield3[0], 1, outfield3[0], 1, outfield3[0], 1); - - // add fields omega = 0.5*(W_x^2 + W_y^2 + W_z^2) - - NekDouble fac = 0.5; - Vmath::Vadd(npoints, &outfield1[0][0], 1, &outfield2[0][0], 1, &omega[0][0], 1); - Vmath::Vadd(npoints, &omega[0][0], 1, &outfield3[0][0], 1, &omega[0][0], 1); - - for (int k = 0; k < addfields; ++k) - { - Vmath::Smul(npoints, fac, &omega[k][0], 1, &omega[k][0], 1); - } - - Vmath::Zero(npoints, &outfield1[0][0], 1); - Vmath::Zero(npoints, &outfield2[0][0], 1); - Vmath::Zero(npoints, &outfield3[0][0], 1); - - Vmath::Vmul(npoints, grad[0 * nfields + 0], 1, grad[0 * nfields + 0], 1, - outfield1[0], 1); - Vmath::Vmul(npoints, grad[1 * nfields + 1], 1, grad[1 * nfields + 1], 1, - outfield2[0], 1); - Vmath::Vmul(npoints, grad[2 * nfields + 2], 1, grad[2 * nfields + 2], 1, - outfield3[0], 1); - - Vmath::Vadd(npoints, &outfield1[0][0], 1, &outfield2[0][0], 1, &S[0][0], 1); - Vmath::Vadd(npoints, &S[0][0], 1, &outfield3[0][0], 1, &S[0][0], 1); - - // W_y + V_z - Vmath::Vadd(npoints, grad[2 * nfields + 1], 1, grad[1 * nfields + 2], 1, - outfield1[0], 1); - Vmath::Vmul(npoints, &outfield1[0][0], 1, &outfield1[0][0], 1, - &outfield1[0][0], 1); - - // U_z + W_x - Vmath::Vadd(npoints, grad[0 * nfields + 2], 1, grad[2 * nfields + 0], 1, - outfield2[0], 1); - Vmath::Vmul(npoints, &outfield2[0][0], 1, &outfield2[0][0], 1, - &outfield2[0][0], 1); - - // V_x + U_y - Vmath::Vadd(npoints, grad[1 * nfields + 0], 1, grad[0 * nfields + 1], 1, - outfield3[0], 1); - Vmath::Vmul(npoints, &outfield3[0][0], 1, &outfield3[0][0], 1, - &outfield3[0][0], 1); - - Vmath::Vadd(npoints, &outfield1[0][0], 1, &outfield2[0][0], 1, - &outfield2[0][0], 1); - Vmath::Vadd(npoints, &outfield2[0][0], 1, &outfield3[0][0], 1, - &outfield3[0][0], 1); - - for (int k = 0; k < addfields; ++k) - { - Vmath::Smul(npoints, fac, &outfield3[k][0], 1, &outfield3[k][0], 1); - } - - Vmath::Vadd(npoints, &outfield3[0][0], 1, &S[0][0], 1, &S[0][0], 1); - Vmath::Vsub(npoints, omega[0], 1, S[0], 1, outfield[0], 1); - - for (int k = 0; k < addfields; ++k) - { - Vmath::Smul(npoints, fac, &outfield[k][0], 1, &outfield[k][0], 1); - } - - + + for (i = 0; i < addfields; ++i) { m_f->m_exp[nfields + i] = m_f->AppendExpList(); m_f->m_exp[nfields + i]->UpdatePhys() = outfield[i]; m_f->m_exp[nfields + i]->FwdTrans(outfield[i], - m_f->m_exp[nfields + i]->UpdateCoeffs()); + m_f->m_exp[nfields + i]->UpdateCoeffs()); } - - vector outname; - outname.push_back("Q"); - + + vector outname; + outname.push_back("Q"); + std::vector FieldDef - = m_f->m_exp[0]->GetFieldDefinitions(); + = m_f->m_exp[0]->GetFieldDefinitions(); std::vector > FieldData(FieldDef.size()); - + for (j = 0; j < nfields + addfields; ++j) { for (i = 0; i < FieldDef.size(); ++i) @@ -227,16 +259,15 @@ namespace Nektar } else { - FieldDef[i]->m_fields.push_back(m_f->m_fielddef[0]->m_fields[j]); + FieldDef[i]->m_fields.push_back( + m_f->m_fielddef[0]->m_fields[j]); } m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]); } } - + m_f->m_fielddef = FieldDef; m_f->m_data = FieldData; - - } } } -- GitLab