109 int main(
int argc,
char *argv[])
154 if(argc >= 6 || argc <4)
156 fprintf(stderr,
"Usage: ./FldCalcBCs meshfile fieldfile streakfile (optional)alpha\n");
159 cout<<
"argc="<<argc<<endl;
161 if(argc==4){ argc=5;}
163 string meshfile(argv[argc-4]);
186 alp_s = argv[argc-1];
187 alpha = boost::lexical_cast<
double>(alp_s);
193 if( vSession->DefinesSolverInfo(
"INTERFACE")
194 &&vSession->GetSolverInfo(
"INTERFACE")==
"phase" )
198 cout<<
"cr="<<cr<<endl;
204 vSession->MatchSolverInfo(
"Symmetrization",
"True",symm,
true);
207 cout<<
"symmetrization is active"<<endl;
212 string fieldfile(argv[argc-3]);
213 vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
214 vector<vector<NekDouble> > fielddata;
219 std::string solvtype = vSession->GetSolverInfo(
"SOLVERTYPE");
221 nfields = fielddef[0]->m_fields.size();
223 if(solvtype ==
"CoupledLinearisedNS" && nfields==1)
230 if(solvtype ==
"CoupledLinearisedNS" && nfields!=2)
232 SetFields(graphShPt,boundaryConditions,vSession,fields,nfields-1);
236 lastfield = nfields-1;
237 cout<<
"Set pressure: "<<lastfield<<endl;
238 int nplanes = fielddef[0]->m_numModes[2];
241 NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
244 fields[lastfield] = Exp3DH1;
247 else if(solvtype ==
"CoupledLinearisedNS" && nfields==2)
249 cout<<
"Set pressure split"<<endl;
251 int nplanes = fielddef[0]->m_numModes[2];
254 NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
255 cout<<
"set ppp"<<endl;
269 lastfield = nfields-1;
272 fields[lastfield] = Exp3DH1;
277 SetFields(graphShPt,boundaryConditions,vSession,fields,nfields);
285 for(
int i = 0; i < fielddata.size(); ++i)
287 fields[0]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[0]->m_fields[0], fields[0]->UpdateCoeffs());
290 fields[0]->BwdTrans_IterPerExp(fields[0]->GetCoeffs(),fields[0]->UpdatePhys());
291 cout<<
"field:"<<fielddef[0]->m_fields[0]<<endl;
302 for(
int i = 0; i < fielddata.size(); ++i)
304 fields[lastfield]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[i]->m_fields[0], fields[lastfield]->UpdateCoeffs());
306 fields[lastfield]->BwdTrans_IterPerExp(fields[lastfield]->GetCoeffs(),fields[lastfield]->UpdatePhys());
319 for(
int j = 0; j < nfields; ++j)
321 for(
int i = 0; i < fielddata.size(); ++i)
323 fields[j]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[i]->m_fields[j], fields[j]->UpdateCoeffs());
325 fields[j]->BwdTrans_IterPerExp(fields[j]->GetCoeffs(),fields[j]->UpdatePhys());
348 string streakfile(argv[argc-2]);
349 vector<LibUtilities::FieldDefinitionsSharedPtr> streakdef;
350 vector<vector<NekDouble> > streakdata;
358 for(
int i = 0; i < streakdata.size(); ++i)
360 streak->ExtractDataToCoeffs(streakdef[i],streakdata[i],streakdef[i]->m_fields[0], streak->UpdateCoeffs());
362 streak->BwdTrans(streak->GetCoeffs(),streak->UpdatePhys());
369 int nIregions, lastIregion;
376 int nbnd= bndConditions.num_elements();
377 for(
int r=0; r<nbnd; r++)
379 if(bndConditions[r]->GetUserDefined()==SpatialDomains::eCalcBC)
383 Ilayers[nIregions]=r;
388 ASSERTL0(nIregions>0,
"there is any boundary region with the tag USERDEFINEDTYPE=""CalcBC"" specified");
437 bool coarseVWI=
false;
443 int coordim = graphShPt->GetMeshDimension();
447 fields, outfieldx,outfieldy,streak,symm,Refindices, alpha,cr);
454 WriteBcs(var,lastIregion, fieldfile,graphShPt,outfieldx);
456 WriteBcs(var,lastIregion, fieldfile,graphShPt,outfieldy);
462 cout<<
"CalcNonLinearForcing"<<endl;
484 bool DeclareCoeffPhysArrays =
true;
500 enum HomogeneousType HomogeneousType = eNotHomogeneous;
502 if(session->DefinesSolverInfo(
"HOMOGENEOUS"))
504 std::string HomoStr = session->GetSolverInfo(
"HOMOGENEOUS");
507 if((HomoStr ==
"HOMOGENEOUS1D")||(HomoStr ==
"Homogeneous1D")||
508 (HomoStr ==
"1D")||(HomoStr ==
"Homo1D"))
510 HomogeneousType = eHomogeneous1D;
511 npointsZ = session->GetParameter(
"HomModesZ");
512 LhomZ = session->GetParameter(
"LZ");
516 if((HomoStr ==
"HOMOGENEOUS2D")||(HomoStr ==
"Homogeneous2D")||
517 (HomoStr ==
"2D")||(HomoStr ==
"Homo2D"))
519 HomogeneousType = eHomogeneous2D;
520 npointsY = session->GetParameter(
"HomModesY");
521 LhomY = session->GetParameter(
"LY");
522 npointsZ = session->GetParameter(
"HomModesZ");
523 LhomZ = session->GetParameter(
"LZ");
527 if((HomoStr ==
"HOMOGENEOUS3D")||(HomoStr ==
"Homogeneous3D")||
528 (HomoStr ==
"3D")||(HomoStr ==
"Homo3D"))
530 HomogeneousType = eHomogeneous3D;
531 npointsX = session->GetParameter(
"HomModesX");
532 LhomX = session->GetParameter(
"LX");
533 npointsY = session->GetParameter(
"HomModesY");
534 LhomY = session->GetParameter(
"LY");
535 npointsZ = session->GetParameter(
"HomModesZ");
536 LhomZ = session->GetParameter(
"LZ");
540 if(session->DefinesSolverInfo(
"USEFFT"))
547 int expdim = mesh->GetMeshDimension();
556 if(HomogeneousType == eHomogeneous2D)
563 for(i = 0 ; i < nvariables; i++)
566 ::AllocateSharedPtr(session,BkeyY,BkeyZ,LhomY,LhomZ,useFFT,deal,mesh,session->GetVariable(i));
571 for(i = 0 ; i < nvariables; i++)
582 if(HomogeneousType == eHomogeneous1D)
586 for(i = 0 ; i < nvariables; i++)
600 for(i = 1 ; i < nvariables; i++)
610 if(HomogeneousType == eHomogeneous3D)
612 ASSERTL0(
false,
"3D fully periodic problems not implemented yet");
622 for(i = 1 ; i < nvariables; i++)
631 ASSERTL0(
false,
"Expansion dimension not recognised");
643 int npts = Ilayer->GetNpoints();
651 Ilayer->GetCoords(coord_x,coord_y);
658 for(i = 0; i <
npts; ++i)
661 xnew = - coord_x[i] + xmax;
664 for(j = start; j >=0 ; --j)
666 if((coord_x[j]-xnew)*(coord_x[j]-xnew) < tol)
677 for(j = npts-1; j > start; --j)
680 if((coord_x[j]-xnew)*(coord_x[j]-xnew) < tol)
686 ASSERTL0(j != start,
"Failsed to find matching point");
703 ASSERTL0(Iregions.num_elements()==3,
"something wrong with the number of I layers");
705 int Ireg =Iregions[2];
706 cout<<
"layer region="<<Ireg<<endl;
713 if(infields.num_elements()==2)
727 int nq1D = Ilayer->GetPlane(0)->GetTotPoints();
733 Ilayer->GetPlane(0)->GetCoords(x0d,x1d,x2d);
746 int np = streak->GetTotPoints();
750 int totcoeffs = streak->GetNcoeffs();
754 streak->PhysDeriv(streak->GetPhys(), gradx, grady);
755 streak->FwdTrans_IterPerExp(gradx, gradxcoeffs);
756 streak->FwdTrans_IterPerExp(grady, gradycoeffs);
818 regIt= bregions.find(Ireg)->second->begin();
819 regIt !=bregions.find(Ireg)->second->end();
822 Icompreg=regIt->second;
829 int nummodes = Ilayer->GetExp(0)->GetNcoeffs();
835 Nregcoeffs = Icompreg->size()*nummodes;
846 int nqedge =nq1D/Icompreg->size();
865 for(
int f=0; f<2; ++f)
871 for(
int k=0; k< Icompreg->size(); k++)
875 !(segmentGeomreg = boost::dynamic_pointer_cast<
878 {
ASSERTL0(
false,
"dynamic cast to a SegGeom failed"); }
880 int EIDreg = segmentGeomreg->GetEid();
881 offsetregIExp = Ilayer->GetPlane(0)->GetCoeff_Offset(k);
884 ->GetElementsFromEdge(segmentGeomreg);
885 int elmtidreg = ((*elementreg)[0]->m_Element)->GetGlobalID();
886 int dimelmtreg = ((*elementreg)[0]->m_Element)->GetNumEdges();
887 Elmtid[k] = elmtidreg;
890 ((*elementreg)[0]->m_Element))->GetEorient((*elementreg)[0]
893 gmat = pressure->GetPlane(0)->GetExp(elmtidreg)->GetMetricInfo()->GetGmat(
894 pressure->GetPlane(0)->GetExp(elmtidreg)->GetPointsKeys());
895 int nq2D = nqedge*nqedge;
898 map<int, int>EdgeGIDreg;
900 for(cnt=f=0; f<dimelmtreg; f++)
903 ((*elementreg)[0]->m_Element))->GetEid(f);
904 EdgeGIDreg[id]=cnt++;
907 int localidedge = EdgeGIDreg.find(EIDreg)->second;
909 localidedge = (*elementreg)[0]->m_EdgeIndx;
916 pressure->GetExp(elmtidreg)->GetEdgeToElementMap(EdgeGIDreg.find(EIDreg)->second,
917 orientreg,bmapedgereg,signreg);
925 int Reoffsetreg = pressure->GetPlane(0)->GetCoeff_Offset(elmtidreg);
926 int Imoffsetreg = pressure->GetPlane(1)->GetCoeff_Offset(elmtidreg);
927 int stoffsetreg = streak->GetCoeff_Offset(elmtidreg);
940 for(
int d=0; d<bmapedgereg.num_elements(); d++)
943 Recoeffsedgereg[d]=Recoeffs[Reoffsetreg+bmapedgereg[d]];
944 Imcoeffsedgereg[d]=Imcoeffs[Imoffsetreg+bmapedgereg[d]];
945 Recoeffsreg[offsetregIExp +d] = Recoeffsedgereg[d];
946 Imcoeffsreg[offsetregIExp +d] = Imcoeffsedgereg[d];
949 stcoeffsedgereg[d] = stcoeffs[stoffsetreg + bmapedgereg[d]];
950 stgradxcoeffsedge[d] = gradxcoeffs[stoffsetreg + bmapedgereg[d]];
951 stgradycoeffsedge[d] = gradycoeffs[stoffsetreg + bmapedgereg[d]];
952 stcoeffsreg[offsetregIExp +d] = stcoeffsedgereg[d];
953 stgradxcoeffsreg[offsetregIExp +d] = stgradxcoeffsedge[d];
954 stgradycoeffsreg[offsetregIExp +d] = stgradycoeffsedge[d];
969 jacedge = Ilayer->GetPlane(0)->GetExp(k)->GetMetricInfo()->GetJac(Ilayer->GetPlane(0)->GetExp(k)->GetPointsKeys());
975 Ilayer->GetPlane(0)->GetExp(k) );
976 edgestdexp->
BwdTrans(Recoeffsedgereg, Pre_edge);
978 edgestdexp->StdPhysDeriv(Pre_edge,
980 for(
int q=0; q<nqedge; q++)
982 Jac[k*nqedge +q]= jacedge[q];
983 dPrestd[k*nqedge +q] = dPre_edge[q];
994 (Ilayer->GetPlane(0)->GetExp(k) );
997 GetLeftAdjacentElementExp()->GetEdgeNormal(localEid);
1002 int physoffsetregIExp = Ilayer->GetPlane(0)->GetPhys_Offset(k);
1003 for(
int e=0; e< nqedge; e++)
1011 nx[physoffsetregIExp +e] = normals[0][e];
1012 ny[physoffsetregIExp +e] = normals[1][e];
1013 tx[physoffsetregIExp +e] = tangents[0][e];
1014 ty[physoffsetregIExp +e] = tangents[1][e];
1026 deriv = pressure->GetPlane(0)->GetExp(elmtidreg)->GetMetricInfo()->GetDeriv(ptsKeys);
1029 int offsetregphys = Ilayer->GetPlane(0)->GetPhys_Offset(k);
1031 Vmath::Vcopy(nq2D, &(deriv[0][0][0]),1, &(derivelmt[0]),1);
1032 pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
1033 Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dxde1[offsetregphys]),1);
1036 Vmath::Vcopy(nq2D, &(deriv[1][1][0]),1, &(derivelmt[0]),1);
1037 pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
1038 Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dyde2[offsetregphys]),1);
1041 Vmath::Vcopy(nq2D, &(deriv[0][1][0]),1, &(derivelmt[0]),1);
1042 pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
1043 Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dyde1[offsetregphys]),1);
1046 Vmath::Vcopy(nq2D, &(deriv[1][0][0]),1, &(derivelmt[0]),1);
1047 pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
1048 Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dxde2[offsetregphys]),1);
1055 pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat0_2D, gmati);
1056 Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat0[offsetregphys]),1);
1058 pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat1_2D, gmati);
1059 Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat1[offsetregphys]),1);
1061 pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat2_2D, gmati);
1062 Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat2[offsetregphys]),1);
1064 pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat3_2D, gmati);
1065 Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat3[offsetregphys]),1);
1069 Ilayer->GetPlane(0)->BwdTrans(Recoeffsreg, Rephysreg);
1070 Ilayer->GetPlane(1)->BwdTrans(Imcoeffsreg, Imphysreg);
1072 Ilayer->GetPlane(0)->BwdTrans(stcoeffsreg, stphysreg);
1073 Ilayer->GetPlane(0)->BwdTrans_IterPerExp(stgradxcoeffsreg, stphysgradxreg);
1074 Ilayer->GetPlane(0)->BwdTrans_IterPerExp(stgradycoeffsreg, stphysgradyreg);
1095 outfieldx->FwdTrans_IterPerExp(dtx, tcoeffs);
1096 outfieldx->BwdTrans_IterPerExp(tcoeffs, dtx);
1098 outfieldx->FwdTrans_IterPerExp(dty, tcoeffs);
1099 outfieldx->BwdTrans_IterPerExp(tcoeffs, dty);
1103 outfieldx->FwdTrans(f_z, fcoeffs);
1104 outfieldx->BwdTrans(fcoeffs, f_z);
1109 outfieldx->FwdTrans(f_zz, fcoeffs);
1110 outfieldx->BwdTrans(fcoeffs, f_zz);
1117 for(
int t=0; t<nq1D; t++)
1120 curv_unsigned[t]=sqrt(dtx[t]*dtx[t] +dty[t]*dty[t]);
1121 delta[t] = 1+f_z[t]*f_z[t];
1123 curv[t] = f_zz[t]/sqrt(delta[t]*delta[t]*delta[t]);
1131 outfieldx->FwdTrans(curv, curv_coeffs);
1132 outfieldx->BwdTrans(curv_coeffs, curv);
1135 for(
int e=0; e< nq1D; e++)
1137 P2reg[e] = Rephysreg[e]*Rephysreg[e] + Imphysreg[e]*Imphysreg[e];
1157 for(
int h=0; h< np; h++)
1159 P2_2D[h] = pressure->GetPlane(0)->GetPhys()[h]*pressure->GetPlane(0)->GetPhys()[h]
1160 + pressure->GetPlane(1)->GetPhys()[h]*pressure->GetPlane(1)->GetPhys()[h];
1162 NekDouble intP2_2D = pressure->GetPlane(0)->PhysIntegral(P2_2D);
1166 NekDouble tmp = pressure->GetPlane(0)->L2(pressure->GetPlane(0)->GetPhys());
1168 tmp = pressure->GetPlane(1)->L2(pressure->GetPlane(1)->GetPhys());
1173 NekDouble Area = pressure->GetPlane(0)->PhysIntegral(I);
1174 norm2D = sqrt(Area/norm2D);
1177 NekDouble Area1 = pressure->GetPlane(0)->PhysIntegral(I_int);
1178 NekDouble normint2D = sqrt(Area1/intP2_2D);
1179 cout<<
"norm2D="<<norm2D<<
" area="<<Area1<<
" intP2_2D="<<intP2_2D<<
" normint2D="<<normint2D<<endl;
1184 for(
int u=0; u<nq1D; u++)
1186 sqrtlen[u] = sqrt(1+f_z[u]*f_z[u]);
1187 s0d[u] = Ilayer->GetPlane(0)->PhysIntegral(sqrtlen);
1189 NekDouble length = Ilayer->GetPlane(0)->PhysIntegral(sqrtlen);
1190 NekDouble int1D = Ilayer->GetPlane(0)->PhysIntegral(P2reg);
1191 norm = sqrt(length/int1D);
1194 cout<<
"norm1D="<<norm<<
" norm2D/norm1D="<<norm2D/norm<<endl;
1195 cout<<
"scal="<<scal<<endl;
1203 for(
int e=0; e< nq1D; e++)
1205 P2reg_aft[e] = Rephysreg[e]*Rephysreg[e] + Imphysreg[e]*Imphysreg[e];
1226 for(
int w=0; w<nq1D; w++)
1228 dUreg[w] = nx[w]*stphysgradxreg[w] +ny[w]*stphysgradyreg[w];
1243 outfieldx->FwdTrans(Jac, jaccoeffs);
1244 outfieldx->BwdTrans(jaccoeffs, Jacreg);
1255 outfieldx->FwdTrans(dP_re, dP_re_coeffs);
1256 outfieldx->FwdTrans(dP_im, dP_im_coeffs);
1257 outfieldx->BwdTrans(dP_re_coeffs, dP_re);
1258 outfieldx->BwdTrans(dP_im_coeffs, dP_im);
1263 cout<<
"x"<<
" P_re"<<
" dP_re"<<
" streak"<<
" dstreak"<<
" pjump"<<endl;
1264 for(
int s=0; s<nq1D; s++)
1266 dP_square[s] = dP_re[s]*dP_re[s] +dP_im[s]*dP_im[s];
1274 outfieldx->FwdTrans(dP_square, dPsquare_coeffs);
1275 outfieldx->BwdTrans(dPsquare_coeffs, dP_square);
1282 double pow = 1.0/3.0;
1284 for(
int y=0; y<nq1D; y++)
1287 mu53[y] = std::pow ((base*base),pow);
1288 mu53[y] = 1/(base*mu53[y]);
1293 outfieldx->FwdTrans(mu53, mucoeffs);
1294 outfieldx->BwdTrans(mucoeffs, mu53);
1299 outfieldx->FwdTrans(prod, prod_coeffs);
1300 outfieldx->BwdTrans(prod_coeffs, prod);
1306 outfieldx->FwdTrans(d2v, prod_coeffs);
1307 outfieldx->BwdTrans(prod_coeffs, d2v);
1403 rho = session->GetParameter(
"RHO");
1407 alpha = session->GetParameter(
"ALPHA");
1409 cout<<
"alpha="<<alpha<<endl;
1411 ASSERTL0(alpha!=0,
"alpha cannot be 0");
1414 alpha53 = std::pow ((alpha*alpha),pow);
1415 alpha53 = 1/(alpha*alpha53);
1416 for(
int c=0; c<nq1D; c++)
1418 pjump[c] = -n0*alpha53*curv[c]*prod[c] ;
1423 cout<<
"alpha^-5/3="<<alpha53<<endl;
1428 outfieldx->FwdTrans(tx, txcoeffs);
1429 outfieldx->FwdTrans(ty, tycoeffs);
1430 outfieldy->BwdTrans(txcoeffs, tx);
1431 outfieldy->BwdTrans(tycoeffs, ty);
1432 cout<<
"RHO=="<<rho<<endl;
1438 for(
int j=0; j<nq1D; j++)
1442 (outfieldx->UpdatePhys())[j] =
1443 rho*rho*(vjump[j]*tx[j]-pjump[j]);
1444 (outfieldy->UpdatePhys())[j] =
1445 rho*rho*(vjump[j]*ty[j]-pjump[j]);
1463 outfieldx->FwdTrans(outfieldx->GetPhys(), finalcoeffs);
1464 outfieldx->BwdTrans(finalcoeffs, outfieldx->UpdatePhys());
1466 outfieldy->FwdTrans(outfieldy->GetPhys(), finalcoeffs);
1467 outfieldy->BwdTrans(finalcoeffs, outfieldy->UpdatePhys());
1475 cout<<
"symmetrise the jump conditions"<<endl;
1477 for(
int i = 0; i < nq1D; ++i)
1479 tmpx[i] =0.5*(outfieldx->GetPhys()[i] - outfieldx->GetPhys()[Refindices[i]]);
1480 tmpy[i] =0.5*(outfieldy->GetPhys()[i] - outfieldy->GetPhys()[Refindices[i]]);
1556 cout<<
"length layer="<<length<<endl;
1566 for(
int g=0; g<nq1D; g++)
1571 invjac2D = (gmat0[g]*tx[g] +gmat2[g]*ty[g]);
1572 jactest = (invjac2D + (1/Jac[g]))/2.;
1573 jactest = (invjac2D- (1/Jac[g]))/jactest;
1574 basetest = abs(curv[g]);
1597 laytest = (abs(stphysreg[g])-cr);
1600 cout<<setw(14)<<x0d[g]<<
" "<<setw(13)<<x1d[g]<<
1601 " "<<Rephysreg[g]<<
" "<<dP_re[g]<<
" "
1602 <<stphysreg[g]<<
" "<<dUreg[g]<<
" "<<mu53[g]<<
" "<<curv[g]
1611 <<setw(13)<<
" "<<dP_square[g]<<
" "<<setw(13)<<prod[g]
1618 tx[g]<<
" "<<ty[g]<<
" "<<
1624 pjump[g]<<setw(13)<<
" "<<d2v[g]<<
" "
1625 <<outfieldx->GetPhys()[g]<<
" "
1626 <<outfieldy->GetPhys()[g]<<
" "<<Imphysreg[g]<<
" "<<dP_im[g]
1632 ASSERTL0(invjac2D*Jac[g]>0,
" sign jac problem..");
1634 ASSERTL0(abs(jactest)<0.3,
"jac 1D problem..");
1644 ASSERTL0(abs(laytest)<0.002,
"critical layer wrong");
1648 cout<<
"WARNING: crit may be wrong"<<endl;
1655 for(
int m=0; m<nq1D; m++)
1657 if(indexjacwarn[m]!=-1)
1659 cout<<
"warning: point with jacerr>20% index="<<m<<
" x="<<x0d[m]<<
" y="<<x1d[m]<<endl;
1663 for(
int a=0; a<Elmtid.num_elements(); a++)
1665 cout<<
"elmt id="<<Elmtid[a]<<
" edge id="<<Edgeid[a]<<endl;
1693 int npts = waveFields[0]->GetPlane(0)->GetNpoints();
1694 int ncoeffs = waveFields[0]->GetPlane(0)->GetNcoeffs();
1703 int nvel = waveFields.num_elements()-1;
1706 for(
int i=0; i< nvel; i++)
1708 waveVelocities[i] = waveFields[i];
1709 Vmath::Vcopy(npts, waveFields[i]->GetPlane(0)->GetPhys(),1, waveVelocities[i]->GetPlane(0)->UpdatePhys(),1 );
1710 Vmath::Vcopy(npts, waveFields[i]->GetPlane(1)->GetPhys(),1, waveVelocities[i]->GetPlane(1)->UpdatePhys(),1 );
1712 Vmath::Vcopy(ncoeffs, waveFields[i]->GetPlane(0)->GetCoeffs(),1, waveVelocities[i]->GetPlane(0)->UpdateCoeffs(),1 );
1713 Vmath::Vcopy(ncoeffs, waveFields[i]->GetPlane(1)->GetCoeffs(),1, waveVelocities[i]->GetPlane(1)->UpdateCoeffs(),1 );
1717 wavePressure = waveFields[nvel];
1720 Vmath::Vcopy(npts, waveFields[nvel]->GetPlane(0)->GetPhys(),1, wavePressure->GetPlane(0)->UpdatePhys(),1 );
1721 Vmath::Vcopy(npts, waveFields[nvel]->GetPlane(1)->GetPhys(),1, wavePressure->GetPlane(1)->UpdatePhys(),1 );
1723 Vmath::Vcopy(ncoeffs, waveFields[nvel]->GetPlane(0)->GetCoeffs(),1, wavePressure->GetPlane(0)->UpdateCoeffs(),1 );
1724 Vmath::Vcopy(ncoeffs, waveFields[nvel]->GetPlane(1)->GetCoeffs(),1, wavePressure->GetPlane(1)->UpdateCoeffs(),1 );
1729 static int projectfield = -1;
1733 if(projectfield == -1)
1737 for(
int i = 0; i < waveVelocities.num_elements(); ++i)
1739 BndConds = waveVelocities[i]->GetBndConditions();
1740 for(
int j = 0; j < BndConds.num_elements(); ++j)
1748 if(projectfield != -1)
1753 if(projectfield == -1)
1755 cout <<
"using first field to project non-linear forcing which imposes a Dirichlet condition" << endl;
1762 wavePressure->GetPlane(0)->BwdTrans(wavePressure->GetPlane(0)->GetCoeffs(),
1763 wavePressure->GetPlane(0)->UpdatePhys());
1764 wavePressure->GetPlane(1)->BwdTrans(wavePressure->GetPlane(1)->GetCoeffs(),
1765 wavePressure->GetPlane(1)->UpdatePhys());
1769 l2 = wavePressure->GetPlane(0)->L2(wavePressure->GetPlane(0)->GetPhys());
1771 l2 = wavePressure->GetPlane(1)->L2(wavePressure->GetPlane(0)->GetPhys());
1774 NekDouble area = waveVelocities[0]->GetPlane(0)->PhysIntegral(der1);
1775 norm = sqrt(area/norm);
1778 waveVelocities[0]->GetPlane(0)->BwdTrans(waveVelocities[0]->GetPlane(0)->GetCoeffs(),waveVelocities[0]->GetPlane(0)->UpdatePhys());
1781 waveVelocities[0]->GetPlane(1)->BwdTrans(waveVelocities[0]->GetPlane(1)->GetCoeffs(),waveVelocities[0]->GetPlane(1)->UpdatePhys());
1784 waveVelocities[1]->GetPlane(0)->BwdTrans(waveVelocities[1]->GetPlane(0)->GetCoeffs(),waveVelocities[1]->GetPlane(0)->UpdatePhys());
1787 waveVelocities[1]->GetPlane(1)->BwdTrans(waveVelocities[1]->GetPlane(1)->GetCoeffs(),waveVelocities[1]->GetPlane(1)->UpdatePhys());
1796 waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
1803 waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
1807 NekDouble rho = session->GetParameter(
"RHO");
1808 NekDouble Re = session->GetParameter(
"RE");
1809 cout<<
"Re="<<Re<<endl;
1810 NekDouble waveForceMag = rho*rho*std::pow(Re,-1./3.);
1815 for(
int i = 1; i < 2; ++i)
1817 vwiForcing[i] = vwiForcing[i-1] + ncoeffs;
1822 waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,vwiForcing[0]);
1826 waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,vwiForcing[0]);
1829 cout<<
"waveforcemag="<<waveForceMag<<endl;
1830 Vmath::Smul(ncoeffs,-waveForceMag, vwiForcing[0],1, vwiForcing[0],1);
1833 waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
1839 waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
1845 waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,vwiForcing[1]);
1849 waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,vwiForcing[1]);
1852 Vmath::Smul(ncoeffs,- waveForceMag, vwiForcing[1],1, vwiForcing[1],1);
1856 cout<<
"symmetrization"<<endl;
1859 waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(vwiForcing[0],der1);
1860 for(
int i = 0; i <
npts; ++i)
1862 val[i] = 0.5*(der1[i] - der1[Refindices[i]]);
1865 waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, vwiForcing[0]);
1868 waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(vwiForcing[1],der1);
1869 for(
int i = 0; i <
npts; ++i)
1871 val[i] = 0.5*(der1[i] - der1[Refindices[i]]);
1874 waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, vwiForcing[1]);
1883 outfield[0] = vwiForcing[0];
1884 outfield[1] = vwiForcing[1];
1886 string sessionName = fieldfile.substr(0,fieldfile.find_last_of(
"."));
1887 std::string outname = sessionName +
".vwi";
1889 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
1890 = waveVelocities[0]->GetPlane(0)->GetFieldDefinitions();
1891 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
1894 for(
int j=0; j< vwiForcing.num_elements(); ++j)
1897 for(
int i=0; i< FieldDef.size(); i++)
1901 FieldDef[i]->m_fields.push_back(var);
1902 waveVelocities[0]->GetPlane(0)->AppendFieldData(FieldDef[i], FieldData[i], vwiForcing[j]);
1915 int npts = wavefield->GetPlane(0)->GetNpoints();
1923 wavefield->GetPlane(0)->GetCoords(coord_x,coord_y);
1929 for(i = 0; i <
npts; ++i)
1931 xnew = - coord_x[i] + xmax;
1932 ynew = - coord_y[i];
1934 for(j = start; j >=0 ; --j)
1936 if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
1947 for(j = npts-1; j > start; --j)
1950 if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
1956 ASSERTL0(j != start,
"Failsed to find matching point");
1965 string outfile = fieldfile.substr(0,fieldfile.find_last_of(
"."));
1966 outfile +=
"_"+variable+
"_";
1968 sprintf(ibnd,
"%d",region);
1970 string endfile(
".bc");
1973 outregionfield->FwdTrans_IterPerExp(outregionfield->GetPhys(),outregionfield->UpdateCoeffs());
1974 fieldcoeffs[0] = outregionfield->UpdateCoeffs();
1975 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
1976 = outregionfield->GetFieldDefinitions();
1977 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
1982 FieldDef[0]->m_fields.push_back(variable);
1984 outregionfield->AppendFieldData(FieldDef[0], FieldData[0]);
1994 if( variables.num_elements()==2)
1999 else if( variables.num_elements()==3)
2007 ASSERTL0(
false,
" something goes wrong... ");
2011 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
2012 = fields[0]->GetFieldDefinitions();
2014 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
2016 for(
int i = 0; i < fields.num_elements(); ++i)
2018 if (fields[i]->GetPhysState()==
true)
2021 fields[i]->FwdTrans(fields[i]->GetPhys(),fields[i]->UpdateCoeffs());
2024 fieldcoeffs[i] = fields[i]->UpdateCoeffs();
2028 for(
int j = 0; j < fieldcoeffs.num_elements(); ++j)
2031 for(
int i = 0; i < FieldDef.size(); ++i)
2034 FieldDef[i]->m_fields.push_back(variables[j]);
2037 fields[0]->AppendFieldData(FieldDef[i], FieldData[i],fieldcoeffs[j]);
void WriteFld(string outfile, SpatialDomains::MeshGraphSharedPtr &mesh, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
#define ASSERTL0(condition, msg)
boost::shared_ptr< ContField1D > ContField1DSharedPtr
std::vector< PointsKey > PointsKeyVector
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.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
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< ContField2D > ContField2DSharedPtr
Array< OneD, int > GetReflectionIndex2D(MultiRegions::ExpListSharedPtr wavefield)
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.
1D Evenly-spaced points using Lagrange polynomial
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
boost::shared_ptr< SegExp > SegExpSharedPtr
1D Evenly-spaced points using Fourier Fit
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< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
void SetFields(SpatialDomains::MeshGraphSharedPtr &mesh, SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions, LibUtilities::SessionReaderSharedPtr &session, Array< OneD, MultiRegions::ExpListSharedPtr > &Exp, int nvariables)
boost::shared_ptr< SegGeom > SegGeomSharedPtr
void Manipulate(Array< OneD, int > Iregions, int coordim, SpatialDomains::MeshGraphSharedPtr &mesh, SpatialDomains::BoundaryConditions &bcs, Array< OneD, MultiRegions::ExpListSharedPtr > &infields, Array< OneD, MultiRegions::ExpList1DSharedPtr > &outfieldx, Array< OneD, MultiRegions::ExpList1DSharedPtr > &outfieldy, MultiRegions::ExpListSharedPtr &streak)
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
boost::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
Defines a specification for a set of points.
void CalcNonLinearForcing(SpatialDomains::MeshGraphSharedPtr &mesh, LibUtilities::SessionReaderSharedPtr &session, string fieldfile, Array< OneD, MultiRegions::ExpListSharedPtr > &waveFields, MultiRegions::ExpListSharedPtr &streak, Array< OneD, int > Refindices, bool symm)
void WriteBcs(string variable, int region, string fieldfile, SpatialDomains::MeshGraphSharedPtr &mesh, MultiRegions::ContField1DSharedPtr &outregionfield)
Array< OneD, int > GetReflectionIndex(MultiRegions::ExpListSharedPtr Exp, int Ireg)
boost::shared_ptr< GeometryVector > Composite
static const NekDouble kNekUnsetDouble
boost::shared_ptr< ElementEdgeVector > ElementEdgeVectorSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void Extractlayerdata(Array< OneD, int > Iregions, int coordim, SpatialDomains::MeshGraphSharedPtr &mesh, LibUtilities::SessionReaderSharedPtr &session, SpatialDomains::BoundaryConditions &bcs, Array< OneD, MultiRegions::ExpListSharedPtr > &infields, MultiRegions::ContField1DSharedPtr &outfieldx, MultiRegions::ContField1DSharedPtr &outfieldy, MultiRegions::ExpListSharedPtr &streak, bool symm, Array< OneD, int > Refindices, NekDouble alpha, NekDouble cr)
boost::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
static const NekDouble kGeomFactorsTol
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space...
int GetLeftAdjacentElementEdge() const
int main(int argc, char *argv[])
const BoundaryRegionCollection & GetBoundaryRegions(void) const
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
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
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Describes the specification for a Basis.
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.