55 ProcessVorticity::create,
"Computes vorticity field.");
69 cout <<
"ProcessVorticity: Calculating vorticity..." << endl;
73 int expdim =
m_f->m_graph->GetMeshDimension();
74 int spacedim = expdim;
75 if ((
m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
76 (
m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
80 int nfields =
m_f->m_fielddef[0]->m_fields.size();
83 ASSERTL0(
false,
"Error: Vorticity for a 1D problem cannot "
86 int addfields = (spacedim == 2)? 1:3;
88 int npoints =
m_f->m_exp[0]->GetNpoints();
94 m_f->m_session->LoadParameter(
"Strip_Z",nstrips,1);
96 m_f->m_exp.resize(nfields*nstrips);
98 for (i = 0; i < spacedim*spacedim; ++i)
103 for (i = 0; i < addfields; ++i)
109 for(
int i = 0; i<spacedim; i++)
114 vector<MultiRegions::ExpListSharedPtr> Exp(nstrips*addfields);
120 for(s = 0; s < nstrips; ++s)
125 if (
m_f->m_fieldMetaDataMap.count(
"MappingCartesianVel"))
127 if(
m_f->m_fieldMetaDataMap[
"MappingCartesianVel"] ==
"False")
130 for (
int i =0; i<spacedim; ++i )
133 if (
m_f->m_exp[0]->GetWaveSpace())
135 m_f->m_exp[0]->HomogeneousBwdTrans(
136 m_f->m_exp[s*nfields+i]->GetPhys(),
147 mapping->ContravarToCartesian(vel, vel);
149 if (
m_f->m_exp[0]->GetWaveSpace())
151 for (
int i =0; i<spacedim; ++i )
153 m_f->m_exp[0]->HomogeneousFwdTrans(vel[i], vel[i]);
159 for (
int i =0; i<spacedim; ++i )
169 for (
int i =0; i<spacedim; ++i )
180 for (i = 0; i < spacedim; ++i)
182 m_f->m_exp[s*nfields+i]->PhysDeriv(vel[i],
185 mapping->CovarToCartesian(tmp, tmp);
186 for(
int j = 0; j<spacedim; j++)
188 Vmath::Vcopy(npoints, tmp[j], 1, grad[i*spacedim+j], 1 );
193 grad[0*spacedim+1], 1,
198 for (i = 0; i < spacedim; ++i)
200 m_f->m_exp[s*nfields+i]->PhysDeriv(vel[i],
204 mapping->CovarToCartesian(tmp, tmp);
205 for(
int j = 0; j<spacedim; j++)
207 Vmath::Vcopy(npoints, tmp[j], 1, grad[i*spacedim+j], 1 );
212 Vmath::Vsub(npoints, grad[2*spacedim+1], 1, grad[1*spacedim+2], 1,
215 Vmath::Vsub(npoints, grad[0*spacedim+2], 1, grad[2*spacedim+0], 1,
218 Vmath::Vsub(npoints, grad[1*spacedim+0], 1, grad[0*spacedim+1], 1,
222 for (i = 0; i < addfields; ++i)
224 int n = s*addfields + i;
225 Exp[n] =
m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir);
226 Exp[n]->UpdatePhys() = outfield[i];
227 Exp[n]->FwdTrans_IterPerExp(outfield[i],
228 Exp[n]->UpdateCoeffs());
233 for(s = 0; s < nstrips; ++s)
235 for(i = 0; i < addfields; ++i)
237 it =
m_f->m_exp.begin()+s*(nfields+addfields)+nfields+i;
238 m_f->m_exp.insert(it, Exp[s*addfields+i]);
242 vector<string > outname;
245 outname.push_back(
"W_z");
249 outname.push_back(
"W_x");
250 outname.push_back(
"W_y");
251 outname.push_back(
"W_z");
254 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
255 =
m_f->m_exp[0]->GetFieldDefinitions();
256 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
258 for(s = 0; s < nstrips; ++s)
260 for (j = 0; j < nfields + addfields; ++j)
262 for (i = 0; i < FieldDef.size()/nstrips; ++i)
264 int n = s * FieldDef.size()/nstrips + i;
268 FieldDef[n]->m_fields.push_back(outname[j-nfields]);
272 FieldDef[n]->m_fields.push_back(
m_f->m_fielddef[0]->m_fields[j]);
274 m_f->m_exp[s*(nfields + addfields)+j]->AppendFieldData(FieldDef[n], FieldData[n]);
279 m_f->m_fielddef = FieldDef;
280 m_f->m_data = FieldData;
#define ASSERTL0(condition, msg)
pair< ModuleType, string > ModuleKey
FieldSharedPtr m_f
Field object.
static GlobalMapping::MappingSharedPtr GetMapping(FieldSharedPtr f)
virtual ~ProcessVorticity()
boost::shared_ptr< Field > FieldSharedPtr
GLOBAL_MAPPING_EXPORT typedef boost::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.