44 std::string ForcingMovingBody::className =
46 "MovingBody", ForcingMovingBody::create,
"Moving Body Forcing");
48 ForcingMovingBody::ForcingMovingBody(
50 const std::weak_ptr<SolverUtils::EquationSystem> &pEquation)
57 const unsigned int &pNumForcingFields,
const TiXmlElement *pForce)
61 "Moving body implemented just for 3D Homogenous 1D expansions.");
94 int phystot = pFields[0]->GetTotPoints();
95 for (
int i = 0; i <
m_zta.size(); i++)
123 int physTot = pFields[0]->GetTotPoints();
126 for (
int i = 0; i < 3; i++)
132 pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
142 m_mapping->UpdateMapping(time, coords, coordsVel);
152 ASSERTL0(
false,
"Motion loading from file needs specific "
153 "implementation: Work in Progress!");
169 int physTot = pFields[0]->GetTotPoints();
172 for (
int i = 0; i < 3; i++)
178 pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
191 for (
int var = 0; var < 3; var++)
193 for (
int plane = 0; plane <
m_np; plane++)
195 int n = pFields[0]->GetPlane(plane)->GetTotPoints();
196 int offset = plane * n;
197 int Offset = var *
m_np + plane;
206 ASSERTL0(
false,
"Unrecogized vibration type for cable's dynamic model");
210 int colrank = vcomm->GetColumnComm()->GetRank();
230 int colrank = vcomm->GetColumnComm()->GetRank();
231 int nproc = vcomm->GetColumnComm()->GetSize();
234 m_session->MatchSolverInfo(
"HomoStrip",
"True", homostrip,
false);
240 npts =
m_session->GetParameter(
"HomModesZ");
244 m_session->LoadParameter(
"HomStructModesZ", npts);
245 m_session->LoadParameter(
"Strip_Z", nstrips);
269 fces[0][0] = Hydroforces[0];
270 fces[1][0] = Hydroforces[
m_np];
276 for (
int i = 1; i < nproc; ++i)
278 vcomm->GetColumnComm()->Recv(i, tmp);
279 for (
int n = 0; n <
m_np; n++)
281 for (
int j = 0; j < 2; ++j)
283 fces[j][i *
m_np + n] = tmp[j *
m_np + n];
293 for (
int i = 1; i < nstrips; ++i)
295 vcomm->GetColumnComm()->Recv(i, tmp);
297 for (
int j = 0; j < 2; ++j)
308 vcomm->GetColumnComm()->Send(0, Hydroforces);
312 for (
int i = 1; i < nstrips; ++i)
317 tmp[0] = Hydroforces[0];
318 tmp[1] = Hydroforces[
m_np];
319 vcomm->GetColumnComm()->Send(0, tmp);
330 m_session->MatchSolverInfo(
"FictitiousMassMethod",
"True", fictmass,
335 m_session->LoadParameter(
"FictMass", fictrho);
336 m_session->LoadParameter(
"FictDamp", fictdamp);
338 static NekDouble Betaq_Coeffs[2][2] = {{1.0, 0.0}, {2.0, -1.0}};
343 int nlevels =
m_fV[0].size();
345 for (
int i = 0; i <
m_motion.size(); ++i)
351 int Aoffset = 2 * npts;
357 Vmath::Smul(npts, Betaq_Coeffs[nint - 1][nint - 1],
358 m_fV[i][nint - 1], 1,
m_fV[i][nlevels - 1], 1);
359 Vmath::Smul(npts, Betaq_Coeffs[nint - 1][nint - 1],
360 m_fA[i][nint - 1], 1,
m_fA[i][nlevels - 1], 1);
362 for (
int n = 0; n < nint - 1; ++n)
365 m_fV[i][nlevels - 1], 1,
m_fV[i][nlevels - 1],
368 m_fA[i][nlevels - 1], 1,
m_fA[i][nlevels - 1],
384 for (
int n = 0, cn = 1; n <
m_vdim; n++, cn--)
398 vcomm->GetColumnComm()->Recv(0, tmp);
403 for (
int i = 1; i < nproc; ++i)
405 for (
int j = 0; j < 2; j++)
407 for (
int k = 0; k < 2; k++)
409 for (
int n = 0; n <
m_np; n++)
416 vcomm->GetColumnComm()->Send(i, tmp);
419 for (
int j = 0; j < 2; j++)
421 for (
int k = 0; k < 2; k++)
423 for (
int n = 0; n <
m_np; n++)
439 for (
int j = 1; j < nproc / nstrips; j++)
441 if (colrank == j * nstrips)
443 vcomm->GetColumnComm()->Recv(0, tmp);
445 for (
int plane = 0; plane <
m_np; plane++)
447 for (
int var = 0; var < 2; var++)
449 for (
int k = 0; k < 2; k++)
451 Motvars[var * 2 *
m_np + k *
m_np + plane] =
459 for (
int i = 1; i < nstrips; i++)
461 for (
int j = 0; j < nproc / nstrips; j++)
463 if (colrank == i + j * nstrips)
465 vcomm->GetColumnComm()->Recv(0, tmp);
467 for (
int plane = 0; plane <
m_np; plane++)
469 for (
int var = 0; var < 2; var++)
471 for (
int k = 0; k < 2; k++)
473 Motvars[var * 2 *
m_np + k *
m_np + plane] =
484 for (
int j = 0; j < 2; ++j)
486 for (
int k = 0; k < 2; ++k)
492 for (
int j = 1; j < nproc / nstrips; j++)
494 vcomm->GetColumnComm()->Send(j * nstrips, tmp);
497 for (
int plane = 0; plane <
m_np; plane++)
499 for (
int j = 0; j < 2; j++)
501 for (
int k = 0; k < 2; ++k)
503 Motvars[j * 2 *
m_np + k *
m_np + plane] =
509 for (
int i = 1; i < nstrips; ++i)
511 for (
int j = 0; j < 2; ++j)
513 for (
int k = 0; k < 2; ++k)
519 for (
int j = 0; j < nproc / nstrips; j++)
521 vcomm->GetColumnComm()->Send(i + j * nstrips, tmp);
528 for (
int var = 0; var < 2; var++)
530 for (
int plane = 0; plane <
m_np; plane++)
532 int n = pFields[0]->GetPlane(plane)->GetTotPoints();
536 int offset = plane * n;
537 int xoffset = var *
m_np + plane;
538 int yoffset = 2 *
m_np + xoffset;
553 std::string supptype =
m_session->GetSolverInfo(
"SupportType");
555 int npts = HydroForces.size();
560 for (
int i = 0; i < 4; i++)
567 for (
int i = 0; i < 3; i++)
569 Vmath::Vcopy(npts, BodyMotions + i * npts, 1, fft_i[i + 1], 1);
573 if (boost::iequals(supptype,
"Free-Free"))
575 for (
int j = 0; j < 4; ++j)
577 m_FFT->FFTFwdTrans(fft_i[j], fft_o[j]);
580 else if (boost::iequals(supptype,
"Pinned-Pinned"))
583 int N = fft_i[0].size();
585 for (
int var = 0; var < 4; var++)
587 for (
int i = 0; i < N; i++)
591 for (
int k = 0; k < N; k++)
594 fft_i[var][k] * sin(M_PI / (N) * (k + 1 / 2) * (i + 1));
601 ASSERTL0(
false,
"Unrecognized support type for cable's motion");
605 for (
int i = 0; i < npts; ++i)
614 for (
int var = 0; var < 3; var++)
616 tmp0[var] = fft_o[var + 1][i];
619 tmp2[0] = fft_o[0][i];
622 nrows, &tmp0[0], 1, 0.0, &tmp1[0], 1);
624 &(
m_CoeffMat_A[i]->GetPtr())[0], nrows, &tmp2[0], 1, 1.0,
627 for (
int var = 0; var < 3; var++)
629 fft_i[var][i] = tmp1[var];
635 if (boost::iequals(supptype,
"Free-Free"))
637 for (
int var = 0; var < 3; var++)
639 m_FFT->FFTBwdTrans(fft_i[var], fft_o[var]);
642 else if (boost::iequals(supptype,
"Pinned-Pinned"))
645 int N = fft_i[0].size();
647 for (
int var = 0; var < 3; var++)
649 for (
int i = 0; i < N; i++)
653 for (
int k = 0; k < N; k++)
655 fft_o[var][i] += fft_i[var][k] *
656 sin(M_PI / (N) * (k + 1) * (i + 1 / 2)) *
664 ASSERTL0(
false,
"Unrecognized support type for cable's motion");
667 for (
int i = 0; i < 3; i++)
670 Vmath::Vcopy(npts, fft_o[i], 1, tmp = BodyMotions + i * npts, 1);
686 int colrank = vcomm->GetColumnComm()->GetRank();
687 int nproc = vcomm->GetColumnComm()->GetSize();
692 m_session->MatchSolverInfo(
"HomoStrip",
"True", homostrip,
false);
696 npts =
m_session->GetParameter(
"HomModesZ");
700 m_session->LoadParameter(
"HomStructModesZ", npts);
708 ZIDs = pFields[0]->GetZIDs();
711 std::string vibtype =
m_session->GetSolverInfo(
"VibrationType");
713 if (boost::iequals(vibtype,
"Constrained"))
717 else if (boost::iequals(vibtype,
"Free"))
721 else if (boost::iequals(vibtype,
"Forced"))
730 int nplanes =
m_session->GetParameter(
"HomModesZ");
739 m_session->LoadParameter(
"DistStrip", DistStrip);
740 m_session->LoadParameter(
"Strip_Z", nstrips);
741 m_lhom = nstrips * DistStrip;
753 m_session->MatchSolverInfo(
"FictitiousMassMethod",
"True", fictmass,
false);
757 m_session->LoadParameter(
"FictMass", fictrho);
758 m_session->LoadParameter(
"FictDamp", fictdamp);
766 for (
int i = 0; i <
m_motion.size(); ++i)
771 for (
int n = 0; n < 2; ++n)
790 std::ifstream inputStream;
795 m_session->LoadParameter(
"HomStructModesZ", nzpoints);
799 nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
802 if (vcomm->GetRank() == 0)
804 std::string filename =
808 inputStream.open(filename.c_str());
812 for (
int n = 0; n < tmp.size(); n++)
814 inputStream >> tmp[n];
819 for (
int n = 0; n < nzpoints; n++)
821 inputStream >> setprecision(6) >> time;
822 inputStream >> setprecision(6) >> z_cds;
824 inputStream >> setprecision(8) >>
826 inputStream >> setprecision(8) >>
829 inputStream >> setprecision(8) >>
831 inputStream >> setprecision(8) >>
843 int nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
846 pFields[0]->GetHomogeneousBasis()->GetZ();
854 for (
int plane = 0; plane <
m_np; plane++)
856 x2[plane] = z_coords[ZIDs[plane]];
866 ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
867 ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
869 int offset = j * npts;
875 for (
int i = 1; i < nproc; ++i)
877 vcomm->GetColumnComm()->Recv(i, tmp0);
878 vcomm->GetColumnComm()->Recv(i, tmp1);
889 vcomm->GetColumnComm()->Send(0, tmp0);
890 vcomm->GetColumnComm()->Send(0, tmp1);
898 m_session->LoadParameter(
"Strip_Z", nstrips);
902 "Fourier transformation of cable motion is currently "
903 "implemented only for FFTW module.");
906 m_session->LoadParameter(
"DistStrip", DistStrip);
914 for (
int plane = 0; plane < npts; plane++)
916 x2[plane] = plane * DistStrip;
921 ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
922 ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
924 int offset = j * npts;
946 m_session->MatchSolverInfo(
"HomoStrip",
"True", homostrip,
false);
950 nplanes =
m_session->GetParameter(
"HomModesZ");
954 m_session->LoadParameter(
"Strip_Z", nplanes);
967 m_session->LoadParameter(
"StructStiff", structstiff, 0.0);
968 m_session->LoadParameter(
"CableTension", cabletension, 0.0);
969 m_session->LoadParameter(
"BendingStiff", bendingstiff, 0.0);
980 std::string supptype =
m_session->GetSolverInfo(
"SupportType");
982 for (
int plane = 0; plane < nplanes; plane++)
994 if (boost::iequals(supptype,
"Free-Free"))
999 else if (boost::iequals(supptype,
"Pinned-Pinned"))
1006 ASSERTL0(
false,
"Unrecognized support type for cable's motion");
1012 tmp7 = tmp2 + tmp4 * tmp6 + tmp5 * tmp7;
1047 int nlevels = input.size();
1050 tmp = input[nlevels - 1];
1052 for (
int n = nlevels - 1; n > 0; --n)
1054 input[n] = input[n - 1];
1079 const TiXmlElement *funcNameElmt_D =
1080 pForce->FirstChildElement(
"DISPLACEMENTS");
1082 "MOVINGBODYFORCE tag has to specify a function name which "
1083 "prescribes the body displacement as d(z,t).");
1087 "Function '" +
m_funcName[0] +
"' not defined.");
1090 const TiXmlElement *funcNameElmt_V =
1091 pForce->FirstChildElement(
"VELOCITIES");
1093 "MOVINGBODYFORCE tag has to specify a function name which "
1094 "prescribes the body velocity of movement as v(z,t).");
1098 "Function '" +
m_funcName[1] +
"' not defined.");
1101 const TiXmlElement *funcNameElmt_A =
1102 pForce->FirstChildElement(
"ACCELERATIONS");
1104 "MOVINGBODYFORCE tag has to specify a function name which "
1105 "prescribes the body acceleration as a(z,t).");
1109 "Function '" +
m_funcName[2] +
"' not defined.");
1125 ASSERTL0(
false,
"The displacements in x must be specified via an "
1126 "equation <E> or a file <F>");
1141 ASSERTL0(
false,
"The displacements in y must be specified via an "
1142 "equation <E> or a file <F>");
1157 ASSERTL0(
false,
"The velocities in x must be specified via an equation "
1158 "<E> or a file <F>");
1173 ASSERTL0(
false,
"The velocities in y must be specified via an equation "
1174 "<E> or a file <F>");
1189 ASSERTL0(
false,
"The accelerations in x must be specified via an "
1190 "equation <E> or a file <F>");
1205 ASSERTL0(
false,
"The accelerations in y must be specified via an "
1206 "equation <E> or a file <F>");
1216 const TiXmlElement *pForce)
1220 std::string typeStr = pForce->Attribute(
"TYPE");
1221 std::map<std::string, std::string> vParams;
1223 const TiXmlElement *param = pForce->FirstChildElement(
"PARAM");
1227 "Missing attribute 'NAME' for parameter in filter " + typeStr +
1229 std::string nameStr = param->Attribute(
"NAME");
1231 ASSERTL0(param->GetText(),
"Empty value string for param.");
1232 std::string valueStr = param->GetText();
1234 vParams[nameStr] = valueStr;
1236 param = param->NextSiblingElement(
"PARAM");
1243 pSession,
m_equ, vParams);
#define ASSERTL0(condition, msg)
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)
GlobalMapping::MappingSharedPtr m_mapping
Array< OneD, std::string > m_funcName
[0] is displacements, [1] is velocities, [2] is accelerations
Array< OneD, Array< OneD, NekDouble > > m_zta
Store the derivatives of motion variables in x-direction.
NekDouble m_lhom
length ratio of the cable
NekDouble m_kinvis
fluid viscous
Array< OneD, DNekMatSharedPtr > m_CoeffMat_A
matrices in Newmart-beta method
Array< OneD, Array< OneD, NekDouble > > m_eta
Store the derivatives of motion variables in y-direction.
void CheckIsFromFile(const TiXmlElement *pForce)
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fV
fictitious velocity storage
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fA
fictitious acceleration storage
void SetDynEqCoeffMatrix(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
void InitialiseCableModel(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
NekDouble m_structrho
mass of the cable per unit length
NekDouble m_timestep
time step
virtual void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
int m_np
number of planes per processors
int m_vdim
vibration dimension
Array< OneD, bool > m_IsFromFile
do determine if the the body motion come from an extern file
FilterMovingBodySharedPtr m_MovBodyfilter
Array< OneD, Array< OneD, NekDouble > > m_MotionVars
storage for the cable's motion(x,y) variables
Array< OneD, DNekMatSharedPtr > m_CoeffMat_B
matrices in Newmart-beta method
void EvaluateStructDynModel(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Hydroforces, NekDouble time)
int m_movingBodyCalls
number of times the movbody have been called
void RollOver(Array< OneD, Array< OneD, NekDouble >> &input)
void Newmark_betaSolver(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &FcePhysinArray, Array< OneD, NekDouble > &MotPhysinArray)
void InitialiseFilter(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
Array< OneD, std::string > m_motion
motion direction: [0] is 'x' and [1] is 'y'
LibUtilities::NektarFFTSharedPtr m_FFT
NekDouble m_structdamp
damping ratio of the cable
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.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Defines a forcing term to be explicitly applied.
const std::weak_ptr< EquationSystem > m_equ
Weak pointer to equation system using this forcing.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, std::string pName, bool pCache=false)
Get a SessionFunction by name.
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
NektarFFTFactory & GetNektarFFTFactory()
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
@ beta
Gauss Radau pinned at x=-1,.
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
@ eFunctionTypeExpression
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
The above copyright notice and this permission notice shall be included.
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 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 Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.