49 const std::weak_ptr<SolverUtils::EquationSystem> &pEquation)
56 [[maybe_unused]]
const unsigned int &pNumForcingFields,
57 const TiXmlElement *pForce)
61 "Moving body implemented just for 3D Homogenous 1D expansions.");
94 size_t phystot = pFields[0]->GetTotPoints();
95 for (
size_t i = 0; i <
m_zta.size(); i++)
124 size_t physTot = pFields[0]->GetTotPoints();
127 for (
size_t i = 0; i < 3; i++)
133 pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
143 m_mapping->UpdateMapping(time, coords, coordsVel);
149 for (
size_t j = 0; j <
m_funcName.size(); j++)
153 ASSERTL0(
false,
"Motion loading from file needs specific "
154 "implementation: Work in Progress!");
170 size_t physTot = pFields[0]->GetTotPoints();
173 for (
size_t i = 0; i < 3; i++)
179 pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
192 for (
size_t var = 0; var < 3; var++)
194 for (
size_t plane = 0; plane <
m_np; plane++)
196 size_t n = pFields[0]->GetPlane(plane)->GetTotPoints();
197 size_t offset = plane * n;
198 size_t Offset = var *
m_np + plane;
207 ASSERTL0(
false,
"Unrecogized vibration type for cable's dynamic model");
211 size_t colrank = vcomm->GetColumnComm()->GetRank();
231 size_t colrank = vcomm->GetColumnComm()->GetRank();
232 size_t nproc = vcomm->GetColumnComm()->GetSize();
235 m_session->MatchSolverInfo(
"HomoStrip",
"True", homostrip,
false);
238 size_t npts, nstrips;
241 npts =
m_session->GetParameter(
"HomModesZ");
245 m_session->LoadParameter(
"HomStructModesZ", npts);
246 m_session->LoadParameter(
"Strip_Z", nstrips);
270 fces[0][0] = Hydroforces[0];
271 fces[1][0] = Hydroforces[
m_np];
277 for (
size_t i = 1; i < nproc; ++i)
279 vcomm->GetColumnComm()->Recv(i, tmp);
280 for (
size_t n = 0; n <
m_np; n++)
282 for (
size_t j = 0; j < 2; ++j)
284 fces[j][i *
m_np + n] = tmp[j *
m_np + n];
294 for (
size_t i = 1; i < nstrips; ++i)
296 vcomm->GetColumnComm()->Recv(i, tmp);
298 for (
size_t j = 0; j < 2; ++j)
309 vcomm->GetColumnComm()->Send(0, Hydroforces);
313 for (
size_t i = 1; i < nstrips; ++i)
318 tmp[0] = Hydroforces[0];
319 tmp[1] = Hydroforces[
m_np];
320 vcomm->GetColumnComm()->Send(0, tmp);
331 m_session->MatchSolverInfo(
"FictitiousMassMethod",
"True", fictmass,
336 m_session->LoadParameter(
"FictMass", fictrho);
337 m_session->LoadParameter(
"FictDamp", fictdamp);
339 static NekDouble Betaq_Coeffs[2][2] = {{1.0, 0.0}, {2.0, -1.0}};
344 size_t nlevels =
m_fV[0].size();
346 for (
size_t i = 0; i <
m_motion.size(); ++i)
351 size_t Voffset = npts;
352 size_t Aoffset = 2 * npts;
358 Vmath::Smul(npts, Betaq_Coeffs[nint - 1][nint - 1],
359 m_fV[i][nint - 1], 1,
m_fV[i][nlevels - 1], 1);
360 Vmath::Smul(npts, Betaq_Coeffs[nint - 1][nint - 1],
361 m_fA[i][nint - 1], 1,
m_fA[i][nlevels - 1], 1);
363 for (
size_t n = 0; n < nint - 1; ++n)
366 m_fV[i][nlevels - 1], 1,
m_fV[i][nlevels - 1],
369 m_fA[i][nlevels - 1], 1,
m_fA[i][nlevels - 1],
385 for (
size_t n = 0, cn = 1; n <
m_vdim; n++, cn--)
399 vcomm->GetColumnComm()->Recv(0, tmp);
404 for (
size_t i = 1; i < nproc; ++i)
406 for (
size_t j = 0; j < 2; j++)
408 for (
size_t k = 0; k < 2; k++)
410 for (
size_t n = 0; n <
m_np; n++)
417 vcomm->GetColumnComm()->Send(i, tmp);
420 for (
size_t j = 0; j < 2; j++)
422 for (
size_t k = 0; k < 2; k++)
424 for (
size_t n = 0; n <
m_np; n++)
440 for (
size_t j = 1; j < nproc / nstrips; j++)
442 if (colrank == j * nstrips)
444 vcomm->GetColumnComm()->Recv(0, tmp);
446 for (
size_t plane = 0; plane <
m_np; plane++)
448 for (
size_t var = 0; var < 2; var++)
450 for (
size_t k = 0; k < 2; k++)
452 Motvars[var * 2 *
m_np + k *
m_np + plane] =
460 for (
size_t i = 1; i < nstrips; i++)
462 for (
size_t j = 0; j < nproc / nstrips; j++)
464 if (colrank == i + j * nstrips)
466 vcomm->GetColumnComm()->Recv(0, tmp);
468 for (
size_t plane = 0; plane <
m_np; plane++)
470 for (
size_t var = 0; var < 2; var++)
472 for (
size_t k = 0; k < 2; k++)
474 Motvars[var * 2 *
m_np + k *
m_np + plane] =
485 for (
size_t j = 0; j < 2; ++j)
487 for (
size_t k = 0; k < 2; ++k)
493 for (
size_t j = 1; j < nproc / nstrips; j++)
495 vcomm->GetColumnComm()->Send(j * nstrips, tmp);
498 for (
size_t plane = 0; plane <
m_np; plane++)
500 for (
size_t j = 0; j < 2; j++)
502 for (
size_t k = 0; k < 2; ++k)
504 Motvars[j * 2 *
m_np + k *
m_np + plane] =
510 for (
size_t i = 1; i < nstrips; ++i)
512 for (
size_t j = 0; j < 2; ++j)
514 for (
size_t k = 0; k < 2; ++k)
520 for (
size_t j = 0; j < nproc / nstrips; j++)
522 vcomm->GetColumnComm()->Send(i + j * nstrips, tmp);
529 for (
size_t var = 0; var < 2; var++)
531 for (
size_t plane = 0; plane <
m_np; plane++)
533 size_t n = pFields[0]->GetPlane(plane)->GetTotPoints();
537 size_t offset = plane * n;
538 size_t xoffset = var *
m_np + plane;
539 size_t yoffset = 2 *
m_np + xoffset;
554 std::string supptype =
m_session->GetSolverInfo(
"SupportType");
556 size_t npts = HydroForces.size();
561 for (
size_t i = 0; i < 4; i++)
568 for (
size_t i = 0; i < 3; i++)
570 Vmath::Vcopy(npts, BodyMotions + i * npts, 1, fft_i[i + 1], 1);
574 if (boost::iequals(supptype,
"Free-Free"))
576 for (
size_t j = 0; j < 4; ++j)
578 m_FFT->FFTFwdTrans(fft_i[j], fft_o[j]);
581 else if (boost::iequals(supptype,
"Pinned-Pinned"))
584 size_t N = fft_i[0].size();
586 for (
size_t var = 0; var < 4; var++)
588 for (
size_t i = 0; i < N; i++)
592 for (
size_t k = 0; k < N; k++)
595 fft_i[var][k] * sin(M_PI / (N) * (k + 1 / 2) * (i + 1));
602 ASSERTL0(
false,
"Unrecognized support type for cable's motion");
606 for (
size_t i = 0; i < npts; ++i)
615 for (
size_t var = 0; var < 3; var++)
617 tmp0[var] = fft_o[var + 1][i];
620 tmp2[0] = fft_o[0][i];
623 nrows, &tmp0[0], 1, 0.0, &tmp1[0], 1);
625 &(
m_CoeffMat_A[i]->GetPtr())[0], nrows, &tmp2[0], 1, 1.0,
628 for (
size_t var = 0; var < 3; var++)
630 fft_i[var][i] = tmp1[var];
636 if (boost::iequals(supptype,
"Free-Free"))
638 for (
size_t var = 0; var < 3; var++)
640 m_FFT->FFTBwdTrans(fft_i[var], fft_o[var]);
643 else if (boost::iequals(supptype,
"Pinned-Pinned"))
646 size_t N = fft_i[0].size();
648 for (
size_t var = 0; var < 3; var++)
650 for (
size_t i = 0; i < N; i++)
654 for (
size_t k = 0; k < N; k++)
656 fft_o[var][i] += fft_i[var][k] *
657 sin(M_PI / (N) * (k + 1) * (i + 1 / 2)) *
665 ASSERTL0(
false,
"Unrecognized support type for cable's motion");
668 for (
size_t i = 0; i < 3; i++)
671 Vmath::Vcopy(npts, fft_o[i], 1, tmp = BodyMotions + i * npts, 1);
687 size_t colrank = vcomm->GetColumnComm()->GetRank();
688 size_t nproc = vcomm->GetColumnComm()->GetSize();
693 m_session->MatchSolverInfo(
"HomoStrip",
"True", homostrip,
false);
697 npts =
m_session->GetParameter(
"HomModesZ");
701 m_session->LoadParameter(
"HomStructModesZ", npts);
709 ZIDs = pFields[0]->GetZIDs();
712 std::string vibtype =
m_session->GetSolverInfo(
"VibrationType");
714 if (boost::iequals(vibtype,
"Constrained"))
718 else if (boost::iequals(vibtype,
"Free"))
722 else if (boost::iequals(vibtype,
"Forced"))
731 int nplanes =
m_session->GetParameter(
"HomModesZ");
740 m_session->LoadParameter(
"DistStrip", DistStrip);
741 m_session->LoadParameter(
"Strip_Z", nstrips);
742 m_lhom = nstrips * DistStrip;
754 m_session->MatchSolverInfo(
"FictitiousMassMethod",
"True", fictmass,
false);
758 m_session->LoadParameter(
"FictMass", fictrho);
759 m_session->LoadParameter(
"FictDamp", fictdamp);
767 for (
size_t i = 0; i <
m_motion.size(); ++i)
772 for (
size_t n = 0; n < 2; ++n)
786 for (
size_t j = 0; j <
m_funcName.size(); j++)
791 std::ifstream inputStream;
796 m_session->LoadParameter(
"HomStructModesZ", nzpoints);
800 nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
803 if (vcomm->GetRank() == 0)
805 std::string filename =
809 inputStream.open(filename.c_str());
813 for (
size_t n = 0; n < tmp.size(); n++)
815 inputStream >> tmp[n];
820 for (
size_t n = 0; n < nzpoints; n++)
822 inputStream >> std::setprecision(6) >> time;
823 inputStream >> std::setprecision(6) >> z_cds;
824 inputStream >> std::setprecision(8) >>
m_MotionVars[0][n];
825 inputStream >> std::setprecision(8) >>
827 inputStream >> std::setprecision(8) >>
829 inputStream >> std::setprecision(8) >>
m_MotionVars[1][n];
830 inputStream >> std::setprecision(8) >>
832 inputStream >> std::setprecision(8) >>
845 pFields[0]->GetHomogeneousBasis()->GetNumModes();
848 pFields[0]->GetHomogeneousBasis()->GetZ();
856 for (
size_t plane = 0; plane <
m_np; plane++)
858 x2[plane] = z_coords[ZIDs[plane]];
868 ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
869 ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
871 size_t offset = j * npts;
877 for (
size_t i = 1; i < nproc; ++i)
879 vcomm->GetColumnComm()->Recv(i, tmp0);
880 vcomm->GetColumnComm()->Recv(i, tmp1);
891 vcomm->GetColumnComm()->Send(0, tmp0);
892 vcomm->GetColumnComm()->Send(0, tmp1);
900 m_session->LoadParameter(
"Strip_Z", nstrips);
904 "Fourier transformation of cable motion is currently "
905 "implemented only for FFTW module.");
908 m_session->LoadParameter(
"DistStrip", DistStrip);
916 for (
size_t plane = 0; plane < npts; plane++)
918 x2[plane] = plane * DistStrip;
923 ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
924 ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
926 size_t offset = j * npts;
948 m_session->MatchSolverInfo(
"HomoStrip",
"True", homostrip,
false);
952 nplanes =
m_session->GetParameter(
"HomModesZ");
956 m_session->LoadParameter(
"Strip_Z", nplanes);
969 m_session->LoadParameter(
"StructStiff", structstiff, 0.0);
970 m_session->LoadParameter(
"CableTension", cabletension, 0.0);
971 m_session->LoadParameter(
"BendingStiff", bendingstiff, 0.0);
982 std::string supptype =
m_session->GetSolverInfo(
"SupportType");
984 for (
size_t plane = 0; plane < nplanes; plane++)
996 if (boost::iequals(supptype,
"Free-Free"))
1001 else if (boost::iequals(supptype,
"Pinned-Pinned"))
1008 ASSERTL0(
false,
"Unrecognized support type for cable's motion");
1014 tmp7 = tmp2 + tmp4 * tmp6 + tmp5 * tmp7;
1049 size_t nlevels = input.size();
1052 tmp = input[nlevels - 1];
1054 for (
size_t n = nlevels - 1; n > 0; --n)
1056 input[n] = input[n - 1];
1081 const TiXmlElement *funcNameElmt_D =
1082 pForce->FirstChildElement(
"DISPLACEMENTS");
1084 "MOVINGBODYFORCE tag has to specify a function name which "
1085 "prescribes the body displacement as d(z,t).");
1089 "Function '" +
m_funcName[0] +
"' not defined.");
1092 const TiXmlElement *funcNameElmt_V =
1093 pForce->FirstChildElement(
"VELOCITIES");
1095 "MOVINGBODYFORCE tag has to specify a function name which "
1096 "prescribes the body velocity of movement as v(z,t).");
1100 "Function '" +
m_funcName[1] +
"' not defined.");
1103 const TiXmlElement *funcNameElmt_A =
1104 pForce->FirstChildElement(
"ACCELERATIONS");
1106 "MOVINGBODYFORCE tag has to specify a function name which "
1107 "prescribes the body acceleration as a(z,t).");
1111 "Function '" +
m_funcName[2] +
"' not defined.");
1127 ASSERTL0(
false,
"The displacements in x must be specified via an "
1128 "equation <E> or a file <F>");
1143 ASSERTL0(
false,
"The displacements in y must be specified via an "
1144 "equation <E> or a file <F>");
1159 ASSERTL0(
false,
"The velocities in x must be specified via an equation "
1160 "<E> or a file <F>");
1175 ASSERTL0(
false,
"The velocities in y must be specified via an equation "
1176 "<E> or a file <F>");
1191 ASSERTL0(
false,
"The accelerations in x must be specified via an "
1192 "equation <E> or a file <F>");
1207 ASSERTL0(
false,
"The accelerations in y must be specified via an "
1208 "equation <E> or a file <F>");
1218 const TiXmlElement *pForce)
1222 std::string typeStr = pForce->Attribute(
"TYPE");
1223 std::map<std::string, std::string> vParams;
1225 const TiXmlElement *param = pForce->FirstChildElement(
"PARAM");
1229 "Missing attribute 'NAME' for parameter in filter " + typeStr +
1231 std::string nameStr = param->Attribute(
"NAME");
1233 ASSERTL0(param->GetText(),
"Empty value string for param.");
1234 std::string valueStr = param->GetText();
1236 vParams[nameStr] = valueStr;
1238 param = param->NextSiblingElement(
"PARAM");
1245 pSession,
m_equ.lock(), vParams);
#define ASSERTL0(condition, msg)
GlobalMapping::MappingSharedPtr m_mapping
Array< OneD, std::string > m_funcName
[0] is displacements, [1] is velocities, [2] is accelerations
ForcingMovingBody(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation)
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
size_t m_np
number of planes per processors
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 RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
void InitialiseCableModel(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
NekDouble m_structrho
mass of the cable per unit length
size_t m_vdim
vibration dimension
NekDouble m_timestep
time step
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) override
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
static SolverUtils::ForcingSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
Creates an instance of this class.
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)
size_t m_movingBodyCalls
number of times the movbody have been called
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
static std::string className
Name of the class.
void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce) override
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 = alpha A x plus beta y 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.
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.