Write mesh to output file.
73 if (
m_f->m_comm->TreatAsRankZero())
75 cout <<
"ProcessScalGrad: Calculating scalar gradient..." << endl;
82 string bvalues =
m_config[
"bnd"].as<
string>();
84 if (bvalues.compare(
"All") == 0)
86 Array<OneD, const MultiRegions::ExpListSharedPtr> BndExp =
87 m_f->m_exp[0]->GetBndCondExpansions();
89 for (i = 0; i < BndExp.num_elements(); ++i)
91 m_f->m_bndRegionsToWrite.push_back(i);
97 m_f->m_bndRegionsToWrite),
98 "Failed to interpret range string");
101 int spacedim =
m_f->m_graph->GetSpaceDimension();
102 if ((
m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
103 (
m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
108 int nfields =
m_f->m_fielddef[0]->m_fields.size();
114 ASSERTL0(
false,
"Error: scalar gradient for a 1D problem cannot "
118 int ngrad = spacedim;
119 int n, cnt, elmtid, nq, offset, boundary, nfq;
120 int npoints =
m_f->m_exp[0]->GetNpoints();
122 Array<OneD, NekDouble> scalar;
123 Array<OneD, Array<OneD, NekDouble> > grad(ngrad), fgrad(ngrad),
128 Array<OneD, int> BoundarytoElmtID, BoundarytoTraceID;
129 Array<OneD, Array<OneD, MultiRegions::ExpListSharedPtr> > BndExp(nfields);
131 m_f->m_exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID, BoundarytoTraceID);
133 for (i = 0; i < nfields; i++)
135 var =
m_f->m_fielddef[0]->m_fields[i];
136 stringstream filename;
137 filename << var <<
"_scalar_gradient";
139 m_f->m_fielddef[0]->m_fields[i] = var;
141 BndExp[i] =
m_f->m_exp[i]->GetBndCondExpansions();
142 outfield[i] = Array<OneD, NekDouble>(npoints);
146 for (cnt = n = 0; n < BndExp[0].num_elements(); ++n)
148 bool doneBnd =
false;
150 for (
int b = 0; b <
m_f->m_bndRegionsToWrite.size(); ++b)
152 if (n ==
m_f->m_bndRegionsToWrite[b])
155 for (i = 0; i < BndExp[0][n]->GetExpSize(); ++i, cnt++)
158 elmtid = BoundarytoElmtID[cnt];
159 elmt =
m_f->m_exp[0]->GetExp(elmtid);
160 nq = elmt->GetTotPoints();
161 offset =
m_f->m_exp[0]->GetPhys_Offset(elmtid);
167 for (j = 0; j < ngrad; ++j)
169 grad[j] = Array<OneD, NekDouble>(nq);
178 for (j = 0; j < nfields; j++)
180 outfield[j] = BndExp[j][n]->UpdateCoeffs() +
181 BndExp[j][n]->GetCoeff_Offset(i);
185 bc = boost::dynamic_pointer_cast<
186 StdRegions::StdExpansion2D>(
187 BndExp[0][n]->GetExp(i));
188 nfq = bc->GetTotPoints();
191 boundary = BoundarytoTraceID[cnt];
195 m_metricinfo = bc->GetMetricInfo();
197 const Array<OneD, const Array<OneD, NekDouble> >
198 normals = elmt->GetFaceNormal(boundary);
201 for (j = 0; j < ngrad; ++j)
203 fgrad[j] = Array<OneD, NekDouble>(nfq);
205 Array<OneD, NekDouble> gradnorm(nfq);
207 for (k = 0; k < nfields; k++)
211 scalar =
m_f->m_exp[k]->GetPhys() + offset;
212 elmt->PhysDeriv(scalar, grad[0], grad[1], grad[2]);
214 for (j = 0; j < ngrad; ++j)
216 elmt->GetFacePhysVals(boundary, bc, grad[j],
221 if (m_metricinfo->GetGtype() ==
224 for (j = 0; j < ngrad; j++)
227 1, gradnorm, 1, gradnorm, 1);
232 for (j = 0; j < ngrad; j++)
235 1, gradnorm, 1, gradnorm, 1);
238 bc->FwdTrans(gradnorm, outfield[k]);
244 if (doneBnd ==
false)
246 cnt += BndExp[0][n]->GetExpSize();
250 for (j = 0; j < nfields; ++j)
252 for (
int b = 0; b <
m_f->m_bndRegionsToWrite.size(); ++b)
254 m_f->m_exp[j]->UpdateBndCondExpansion(
m_f->m_bndRegionsToWrite[b]) =
255 BndExp[j][
m_f->m_bndRegionsToWrite[b]];
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
static bool GenerateOrderedVector(const char *const str, std::vector< unsigned int > &vec)
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< StdExpansion2D > StdExpansion2DSharedPtr
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Geometry is curved or has non-constant factors.
FieldSharedPtr m_f
Field object.