45 RegisterCreatorFunction(
"MovingBody",
46 ForcingMovingBody::create,
47 "Moving Body Forcing");
49 ForcingMovingBody::ForcingMovingBody(
57 const unsigned int& pNumForcingFields,
58 const TiXmlElement* pForce)
62 "Moving body implemented just for 3D Homogenous 1D expansions.");
94 int phystot = pFields[0]->GetTotPoints();
95 for(
int i = 0; i <
m_zta.num_elements(); i++)
124 int physTot = pFields[0]->GetTotPoints();
127 for(
int i =0; i<3; i++)
133 pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
143 m_mapping->UpdateMapping(time, coords, coordsVel);
149 for(
int j = 0; j <
m_funcName.num_elements(); j++)
153 ASSERTL0(
false,
"Motion loading from file needs specific "
154 "implementation: Work in Progress!");
170 int physTot = pFields[0]->GetTotPoints();
173 for(
int i =0; i<3; i++)
179 pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
192 for(
int var = 0; var < 3; var++)
194 for(
int plane = 0; plane <
m_np; plane++)
196 int n = pFields[0]->GetPlane(plane)->GetTotPoints();
197 int offset = plane * n;
198 int Offset = var * m_np+plane;
208 "Unrecogized vibration type for cable's dynamic model");
212 int colrank = vcomm->GetColumnComm()->GetRank();
235 int colrank = vcomm->GetColumnComm()->GetRank();
236 int nproc = vcomm->GetColumnComm()->GetSize();
239 m_session->MatchSolverInfo(
"HomoStrip",
"True",homostrip,
false);
245 npts =
m_session->GetParameter(
"HomModesZ");
249 m_session->LoadParameter(
"HomStructModesZ", npts);
250 m_session->LoadParameter(
"Strip_Z", nstrips);
273 fces[0][0] = Hydroforces[0];
274 fces[1][0] = Hydroforces[
m_np];
280 for (
int i = 1; i < nproc; ++i)
282 vcomm->GetColumnComm()->Recv(i, tmp);
283 for(
int n = 0; n <
m_np; n++)
285 for(
int j = 0; j < 2; ++j)
287 fces[j][i*m_np + n] = tmp[j*m_np + n];
297 for(
int i = 1; i < nstrips; ++i)
299 vcomm->GetColumnComm()->Recv(i, tmp);
301 for(
int j = 0 ; j < 2; ++j)
312 vcomm->GetColumnComm()->Send(0, Hydroforces);
316 for(
int i = 1; i < nstrips; ++i)
321 tmp[0] = Hydroforces[0];
322 tmp[1] = Hydroforces[
m_np];
323 vcomm->GetColumnComm()->Send(0, tmp);
334 m_session->MatchSolverInfo(
"FictitiousMassMethod",
"True",
339 m_session->LoadParameter(
"FictMass", fictrho);
340 m_session->LoadParameter(
"FictDamp", fictdamp);
343 {{1.0, 0.0},{2.0, -1.0}};
348 int nlevels =
m_fV[0].num_elements();
350 for(
int i = 0; i <
m_motion.num_elements(); ++i)
356 int Aoffset = 2*
npts;
363 Betaq_Coeffs[nint-1][nint-1],
365 m_fV[i][nlevels-1], 1);
367 Betaq_Coeffs[nint-1][nint-1],
369 m_fA[i][nlevels-1], 1);
371 for(
int n = 0; n < nint-1; ++n)
374 Betaq_Coeffs[nint-1][n],
376 m_fV[i][nlevels-1],1);
378 Betaq_Coeffs[nint-1][n],
380 m_fA[i][nlevels-1],1);
385 fces[i],1,fces[i],1);
387 fces[i],1,fces[i],1);
395 for(
int n = 0, cn = 1; n <
m_vdim; n++, cn--)
409 vcomm->GetColumnComm()->Recv(0, tmp);
414 for (
int i = 1; i < nproc; ++i)
416 for(
int j = 0; j < 2; j++)
418 for(
int k = 0; k < 2; k++)
420 for (
int n = 0; n <
m_np; n++)
422 tmp[j*2*m_np+k*m_np+n] =
m_MotionVars[j][k*npts+i*m_np+n];
426 vcomm->GetColumnComm()->Send(i, tmp);
429 for(
int j = 0; j < 2; j++)
431 for(
int k = 0; k < 2; k++)
433 for(
int n = 0; n <
m_np; n++)
448 for (
int j = 1; j < nproc/nstrips; j++)
450 if(colrank == j*nstrips)
452 vcomm->GetColumnComm()->Recv(0, tmp);
454 for(
int plane = 0; plane <
m_np; plane++)
456 for(
int var = 0; var < 2; var++)
458 for(
int k = 0; k < 2; k++)
460 Motvars[var*2*m_np+k*m_np+plane]= tmp[var*2+k];
467 for(
int i = 1; i < nstrips; i++)
469 for (
int j = 0; j < nproc/nstrips; j++)
471 if(colrank == i+j*nstrips)
473 vcomm->GetColumnComm()->Recv(0, tmp);
475 for(
int plane = 0; plane <
m_np; plane++)
477 for(
int var = 0; var < 2; var++)
479 for(
int k = 0; k < 2; k++)
481 Motvars[var*2*m_np+k*m_np+plane] = tmp[var*2+k];
491 for(
int j = 0; j < 2; ++j)
493 for(
int k = 0; k < 2; ++k)
499 for (
int j = 1; j < nproc/nstrips; j++)
501 vcomm->GetColumnComm()->Send(j*nstrips, tmp);
504 for(
int plane = 0; plane <
m_np; plane++)
506 for(
int j = 0; j < 2; j++)
508 for(
int k = 0; k < 2; ++k)
515 for(
int i = 1; i < nstrips; ++i)
517 for(
int j = 0; j < 2; ++j)
519 for(
int k = 0; k < 2; ++k)
525 for (
int j = 0; j < nproc/nstrips; j++)
527 vcomm->GetColumnComm()->Send(i+j*nstrips, tmp);
534 for(
int var = 0; var < 2; var++)
536 for(
int plane = 0; plane <
m_np; plane++)
538 int n = pFields[0]->GetPlane(plane)->GetTotPoints();
542 int offset = plane * n;
543 int xoffset = var * m_np+plane;
544 int yoffset = 2*m_np + xoffset;
561 std::string supptype =
m_session->GetSolverInfo(
"SupportType");
563 int npts = HydroForces.num_elements();
568 for(
int i = 0; i < 4; i++)
575 for(
int i = 0; i < 3; i++)
577 Vmath::Vcopy(npts, BodyMotions+i*npts, 1, fft_i[i+1], 1);
581 if(boost::iequals(supptype,
"Free-Free"))
583 for(
int j = 0 ; j < 4; ++j)
585 m_FFT->FFTFwdTrans(fft_i[j], fft_o[j]);
588 else if(boost::iequals(supptype,
"Pinned-Pinned"))
591 int N = fft_i[0].num_elements();
593 for(
int var = 0; var < 4; var++)
595 for(
int i = 0; i < N; i++)
599 for(
int k = 0; k < N; k++)
603 sin(M_PI/(N)*(k+1/2)*(i+1));
611 "Unrecognized support type for cable's motion");
615 for(
int i = 0; i <
npts; ++i)
624 for(
int var = 0; var < 3; var++)
626 tmp0[var] = fft_o[var+1][i];
629 tmp2[0] = fft_o[0][i];
631 Blas::Dgemv(
'N', nrows, nrows, 1.0,
640 for(
int var = 0; var < 3; var++)
642 fft_i[var][i] = tmp1[var];
648 if(boost::iequals(supptype,
"Free-Free"))
650 for(
int var = 0; var < 3; var++)
652 m_FFT->FFTBwdTrans(fft_i[var], fft_o[var]);
655 else if(boost::iequals(supptype,
"Pinned-Pinned"))
658 int N = fft_i[0].num_elements();
660 for(
int var = 0; var < 3; var++)
662 for(
int i = 0; i < N; i++)
666 for(
int k = 0; k < N; k++)
670 sin(M_PI/(N)*(k+1)*(i+1/2))*2/N;
678 "Unrecognized support type for cable's motion");
682 for(
int i = 0; i < 3; i++)
685 Vmath::Vcopy(npts, fft_o[i], 1, tmp = BodyMotions+i*npts, 1);
701 int colrank = vcomm->GetColumnComm()->GetRank();
702 int nproc = vcomm->GetColumnComm()->GetSize();
707 m_session->MatchSolverInfo(
"HomoStrip",
"True",homostrip,
false);
711 npts =
m_session->GetParameter(
"HomModesZ");
715 m_session->LoadParameter(
"HomStructModesZ", npts);
723 ZIDs = pFields[0]->GetZIDs();
724 m_np = ZIDs.num_elements();
726 std::string vibtype =
m_session->GetSolverInfo(
"VibrationType");
728 if(boost::iequals(vibtype,
"Constrained"))
732 else if (boost::iequals(vibtype,
"Free"))
736 else if (boost::iequals(vibtype,
"Forced"))
745 int nplanes =
m_session->GetParameter(
"HomModesZ");
755 m_session->LoadParameter(
"DistStrip", DistStrip);
756 m_session->LoadParameter(
"Strip_Z", nstrips);
757 m_lhom = nstrips * DistStrip;
770 m_session->MatchSolverInfo(
"FictitiousMassMethod",
"True",
775 m_session->LoadParameter(
"FictMass", fictrho);
776 m_session->LoadParameter(
"FictDamp", fictdamp);
784 for (
int i = 0; i <
m_motion.num_elements(); ++i)
789 for(
int n = 0; n < 2; ++n)
803 for(
int j = 0; j <
m_funcName.num_elements(); j++)
808 std::ifstream inputStream;
813 m_session->LoadParameter(
"HomStructModesZ", nzpoints);
817 nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
820 if (vcomm->GetRank() == 0)
822 std::string filename =
m_session->GetFunctionFilename(
826 inputStream.open(filename.c_str());
830 for(
int n = 0; n < tmp.num_elements(); n++)
832 inputStream >> tmp[n];
837 for (
int n = 0; n < nzpoints; n++)
839 inputStream >> setprecision(6) >> time;
840 inputStream >> setprecision(6) >> z_cds;
841 inputStream >> setprecision(8) >> m_MotionVars[0][n];
842 inputStream >> setprecision(8) >> m_MotionVars[0][n+nzpoints];
843 inputStream >> setprecision(8) >> m_MotionVars[0][n+2*nzpoints];
844 inputStream >> setprecision(8) >> m_MotionVars[1][n];
845 inputStream >> setprecision(8) >> m_MotionVars[1][n+nzpoints];
846 inputStream >> setprecision(8) >> m_MotionVars[1][n+2*nzpoints];
857 int nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
860 = pFields[0]->GetHomogeneousBasis()->GetZ();
868 for (
int plane = 0; plane <
m_np; plane++)
870 x2[plane] = z_coords[ZIDs[plane]];
880 ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
881 ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
884 Vmath::Vcopy(m_np, tmp0, 1, tmp = m_MotionVars[0]+offset,1);
885 Vmath::Vcopy(m_np, tmp1, 1, tmp = m_MotionVars[1]+offset,1);
889 for (
int i = 1; i < nproc; ++i)
891 vcomm->GetColumnComm()->Recv(i, tmp0);
892 vcomm->GetColumnComm()->Recv(i, tmp1);
893 Vmath::Vcopy(m_np, tmp0, 1, tmp = m_MotionVars[0]+offset+i*m_np,1);
894 Vmath::Vcopy(m_np, tmp1, 1, tmp = m_MotionVars[1]+offset+i*m_np,1);
899 vcomm->GetColumnComm()->Send(0, tmp0);
900 vcomm->GetColumnComm()->Send(0, tmp1);
908 m_session->LoadParameter(
"Strip_Z", nstrips);
911 "Fourier transformation of cable motion is currently "
912 "implemented only for FFTW module.");
915 m_session->LoadParameter(
"DistStrip", DistStrip);
923 for (
int plane = 0; plane <
npts; plane++)
925 x2[plane] = plane*DistStrip;
930 ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
931 ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
934 Vmath::Vcopy(npts, tmp0, 1, tmp = m_MotionVars[0]+offset,1);
935 Vmath::Vcopy(npts, tmp1, 1, tmp = m_MotionVars[1]+offset,1);
954 m_session->MatchSolverInfo(
"HomoStrip",
"True",homostrip,
false);
958 nplanes =
m_session->GetParameter(
"HomModesZ");
962 m_session->LoadParameter(
"Strip_Z", nplanes);
975 m_session->LoadParameter(
"StructStiff", structstiff, 0.0);
976 m_session->LoadParameter(
"CableTension", cabletension, 0.0);
977 m_session->LoadParameter(
"BendingStiff", bendingstiff, 0.0);
988 std::string supptype =
m_session->GetSolverInfo(
"SupportType");
990 for(
int plane = 0; plane < nplanes; plane++)
1002 if (boost::iequals(supptype,
"Free-Free"))
1005 beta = 2.0 * M_PI/
m_lhom;
1007 else if(boost::iequals(supptype,
"Pinned-Pinned"))
1015 "Unrecognized support type for cable's motion");
1021 tmp7 = tmp2 + tmp4 * tmp6 + tmp5 * tmp7;
1032 (*m_CoeffMat_B[plane])(0,0) = 0.0;
1033 (*m_CoeffMat_B[plane])(0,1) = 0.0;
1034 (*m_CoeffMat_B[plane])(0,2) = 0.0;
1035 (*m_CoeffMat_B[plane])(1,0) = 0.0;
1036 (*m_CoeffMat_B[plane])(1,1) = 1.0;
1037 (*m_CoeffMat_B[plane])(1,2) = m_timestep/2.0;
1038 (*m_CoeffMat_B[plane])(2,0) = 1.0;
1040 (*m_CoeffMat_B[plane])(2,2) = tmp1/4.0;
1043 (*m_CoeffMat_B[plane]) =
1057 int nlevels = input.num_elements();
1060 tmp = input[nlevels-1];
1062 for(
int n = nlevels-1; n > 0; --n)
1064 input[n] = input[n-1];
1089 const TiXmlElement* funcNameElmt_D
1090 = pForce->FirstChildElement(
"DISPLACEMENTS");
1092 "MOVINGBODYFORCE tag has to specify a function name which "
1093 "prescribes the body displacement as d(z,t).");
1097 "Function '" +
m_funcName[0] +
"' not defined.");
1100 const TiXmlElement* funcNameElmt_V
1101 = pForce->FirstChildElement(
"VELOCITIES");
1103 "MOVINGBODYFORCE tag has to specify a function name which "
1104 "prescribes the body velocity of movement as v(z,t).");
1108 "Function '" +
m_funcName[1] +
"' not defined.");
1112 const TiXmlElement* funcNameElmt_A
1113 = pForce->FirstChildElement(
"ACCELERATIONS");
1115 "MOVINGBODYFORCE tag has to specify a function name which "
1116 "prescribes the body acceleration as a(z,t).");
1120 "Function '" +
m_funcName[2] +
"' not defined.");
1136 ASSERTL0(
false,
"The displacements in x must be specified via an "
1137 "equation <E> or a file <F>");
1152 ASSERTL0(
false,
"The displacements in y must be specified via an "
1153 "equation <E> or a file <F>");
1168 ASSERTL0(
false,
"The velocities in x must be specified via an equation "
1169 "<E> or a file <F>");
1184 ASSERTL0(
false,
"The velocities in y must be specified via an equation "
1185 "<E> or a file <F>");
1200 ASSERTL0(
false,
"The accelerations in x must be specified via an "
1201 "equation <E> or a file <F>");
1216 ASSERTL0(
false,
"The accelerations in y must be specified via an "
1217 "equation <E> or a file <F>");
1227 const TiXmlElement* pForce)
1231 std::string typeStr = pForce->Attribute(
"TYPE");
1232 std::map<std::string, std::string> vParams;
1234 const TiXmlElement *param = pForce->FirstChildElement(
"PARAM");
1238 "Missing attribute 'NAME' for parameter in filter "
1240 std::string nameStr = param->Attribute(
"NAME");
1242 ASSERTL0(param->GetText(),
"Empty value string for param.");
1243 std::string valueStr = param->GetText();
1245 vParams[nameStr] = valueStr;
1247 param = param->NextSiblingElement(
"PARAM");
Array< OneD, DNekMatSharedPtr > m_CoeffMat_B
matrices in Newmart-beta method
NekDouble m_structrho
mass of the cable per unit length
#define ASSERTL0(condition, msg)
void CheckIsFromFile(const TiXmlElement *pForce)
NekDouble m_structdamp
damping ratio of the cable
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Array< OneD, Array< OneD, NekDouble > > m_MotionVars
storage for the cable's motion(x,y) variables
void SetDynEqCoeffMatrix(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Array< OneD, DNekMatSharedPtr > m_CoeffMat_A
matrices in Newmart-beta method
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
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
virtual void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, MultiRegions::ExpListSharedPtr > pFields, LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))
GlobalMapping::MappingSharedPtr m_mapping
void Newmark_betaSolver(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &FcePhysinArray, Array< OneD, NekDouble > &MotPhysinArray)
NektarFFTFactory & GetNektarFFTFactory()
void InitialiseFilter(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
int m_movingBodyCalls
number of times the movbody have been called
NekDouble m_kinvis
fluid viscous
Array< OneD, std::string > m_motion
motion direction: [0] is 'x' and [1] is 'y'
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
int m_vdim
vibration dimension
LibUtilities::NektarFFTSharedPtr m_FFT
virtual void v_Apply(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
NekDouble m_timestep
time step
boost::shared_ptr< Equation > EquationSharedPtr
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.
void InitialiseCableModel(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Array< OneD, bool > m_IsFromFile
do determine if the the body motion come from an extern file
Array< OneD, Array< OneD, NekDouble > > m_eta
Store the derivatives of motion variables in y-direction.
FilterMovingBodySharedPtr m_MovBodyfilter
NekDouble m_lhom
length ratio of the cable
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fV
fictitious velocity storage
int m_np
number of planes per processors
Array< OneD, std::string > m_funcName
[0] is displacements, [1] is velocities, [2] is accelerations
void EvaluateStructDynModel(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Hydroforces, NekDouble time)
Defines a forcing term to be explicitly applied.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fA
fictitious acceleration storage
void RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Array< OneD, Array< OneD, NekDouble > > m_zta
Store the derivatives of motion variables in x-direction.
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
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.