13 int main(
int argc,
char *argv[])
20 fprintf(stderr,
"Usage: FldAddScalGrad meshfile infld outfld BoundaryID\n");
24 surfID = boost::lexical_cast<
int>(argv[argc - 1]);
26 argv[argc -1] = argv[argc - 2];
33 string meshfile(argv[argc-4]);
39 string fieldfile(argv[argc-3]);
40 vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
41 vector<vector<NekDouble> > fielddata;
47 int expdim = graphShPt->GetMeshDimension();
57 ASSERTL0(
false,
"Expansion dimension not recognised");
62 ASSERTL0(
false,
"Expansion dimension not recognised");
70 vSession->GetVariable(0));
72 m_locToGlobalMap = originalfield->GetLocalToGlobalMap();
74 exp[0] = originalfield;
75 for (i=0; i<expdim; i++)
81 for (i=expdim; i<addfields; i++)
85 vSession->GetVariable(0));
91 ASSERTL0(
false,
"Expansion dimension not recognised");
98 for(j = 0; j < nfields+addfields; ++j)
100 for(
int i = 0; i < fielddata.size(); ++i)
102 exp[j]->ExtractDataToCoeffs(fielddef [i],
104 fielddef [i]->m_fields[0],
105 exp[j]->UpdateCoeffs());
108 exp[j]->BwdTrans(exp[j]->GetCoeffs(),exp[j]->UpdatePhys());
115 int n, cnt, elmtid, nq, offset, nt, boundary, nfq;
116 nt = exp[0]->GetNpoints();
130 exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
133 for (i = 0; i<addfields-expdim; i++)
135 BndExp[i] = exp[i+1+expdim]->GetBndCondExpansions();
139 exp[0]->PhysDeriv(exp[0]->GetPhys(),exp[1]->UpdatePhys(),
140 exp[2]->UpdatePhys(),exp[3]->UpdatePhys());
144 for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
149 for(i = 0; i < BndExp[0][n]->GetExpSize(); ++i, cnt++)
152 elmtid = BoundarytoElmtID[cnt];
153 elmt = exp[0]->GetExp(elmtid);
154 nq = elmt->GetTotPoints();
155 offset = exp[0]->GetPhys_Offset(elmtid);
159 for(j = 0; j < expdim; ++j)
169 for (j = 0; j< addfields-expdim; j++)
171 values[j] = BndExp[j][n]->UpdateCoeffs() + BndExp[0][n]->GetCoeff_Offset(i);
181 boundary = BoundarytoTraceID[cnt];
184 U = exp[0]->GetPhys() + offset;
187 elmt->PhysDeriv(U,grad[0],grad[1],grad[2]);
190 for(j = 0; j < expdim; ++j)
197 for(j = 0; j < expdim; ++j)
199 elmt->GetFacePhysVals(boundary,bc,grad[j],fgrad[j]);
206 = elmt->GetFaceNormal(boundary);
213 normals[1],1,fgrad[1],1,gradnorm,1);
214 Vmath::Vvtvp (nfq,normals[2],1,fgrad[2],1,gradnorm,1,gradnorm,1);
219 normals[1][0],fgrad[1],1,gradnorm,1);
220 Vmath::Svtvp(nfq,normals[2][0],fgrad[2],1,gradnorm,1,gradnorm,1);
223 for(j = 0; j<expdim; j++)
225 bc->FwdTrans(normals[j],values[j]);
230 bc->FwdTrans(gradnorm,values[expdim]);
237 cnt += BndExp[0][n]->GetExpSize();
241 for(j = 0; j< expdim; j++)
243 exp[j+1]->FwdTrans(exp[j+1]->GetPhys(), exp[j+1]->UpdateCoeffs());
246 for(j = 0; j < addfields-expdim; ++j)
248 int ncoeffs = exp[0]->GetNcoeffs();
251 output=exp[j+1+expdim]->UpdateCoeffs();
253 int nGlobal=m_locToGlobalMap->GetNumGlobalCoeffs();
261 for(
int i = 0; i < BndExp[j].num_elements(); ++i)
266 for(
int k = 0; k < (BndExp[j][i])->GetNcoeffs(); ++k)
268 sign = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
269 outarray[map[bndcnt++]] = sign * coeffs[k];
274 bndcnt += BndExp[j][i]->GetNcoeffs();
277 m_locToGlobalMap->GlobalToLocal(outarray,output);
283 string out(argv[argc-2]);
284 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
285 = exp[0]->GetFieldDefinitions();
286 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
288 vector<string > outname;
290 outname.push_back(
"du/dx");
291 outname.push_back(
"du/dy");
292 outname.push_back(
"du/dz");
293 outname.push_back(
"nx");
294 outname.push_back(
"ny");
295 outname.push_back(
"nz");
296 outname.push_back(
"gradient");
299 for(j = 0; j < nfields+addfields; ++j)
301 for(i = 0; i < FieldDef.size(); ++i)
305 FieldDef[i]->m_fields.push_back(outname[j-nfields]);
309 FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
311 exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
#define ASSERTL0(condition, msg)
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > ElementiDs)
Imports an FLD file.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap)
Write a field file in serial only.
boost::shared_ptr< ContField3D > ContField3DSharedPtr
boost::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Geometry is curved or has non-constant factors.
int main(int argc, char *argv[])