diff --git a/examples/perpendicular-flap-multipatch-vertex-gismo.cpp b/examples/perpendicular-flap-multipatch-vertex-gismo.cpp new file mode 100644 index 0000000..215edc2 --- /dev/null +++ b/examples/perpendicular-flap-multipatch-vertex-gismo.cpp @@ -0,0 +1,466 @@ +/** @file perpendicular-flap-multipatch-vertex-gismo.cpp + + @brief Elasticity participant for the PreCICE example "perpendicular-flap" using multi-patches. + + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): J.Li (2023 - ..., TU Delft), H.M. Verhelst (2019 - 2024, TU Delft) +*/ + +#include +#include +#include +#include +#include + +#include +#include + +#ifdef gsStructuralAnalysis_ENABLED +#include +#include +#include +#include +#include +#include +#include +#endif + +using namespace gismo; + +int main(int argc, char *argv[]) +{ + //! [Parse command line] + bool plot = false; + index_t plotmod = 10; + index_t numRefine = 0; + index_t numElevate = 0; + std::string precice_config; + bool nonlinear = false; + int method = 3; // 1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson, 6 RK4 + + gsCmdLine cmd("Flow over heated plate for PreCICE."); + cmd.addInt( "e", "degreeElevation", + "Number of degree elevation steps to perform before solving (0: equalize degree in all directions)", numElevate ); + cmd.addInt( "r", "uniformRefine", "Number of Uniform h-refinement loops", numRefine ); + cmd.addString( "c", "config", "PreCICE config file", precice_config ); + cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot); + cmd.addInt("m","plotmod", "Modulo for plotting, i.e. if plotmod==1, plots every timestep", plotmod); + cmd.addInt("M", "method","1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson, 6: RK4",method); + cmd.addSwitch("nonlinear", "Use nonlinear elasticity", nonlinear); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + //! [Read input file] + GISMO_ASSERT(gsFileManager::fileExists(precice_config),"No precice config file has been defined"); + + // Generate domain + gsMultiPatch<> patches; + // get from XML + + gsKnotVector<> KV1 (0,1,2,2) ; + gsKnotVector<> KV2 (0,1,12,2) ; + + gsTensorBSplineBasis<2> basis(KV1,KV2); + gsMatrix<> coefs = basis.anchors(); + coefs.transposeInPlace(); + coefs.col(0).array() *= 0.1; + coefs.col(0).array() -= 0.05; + coefs.col(1).array() *= 0.5; + + + patches.addPatch(basis.makeGeometry(coefs)); + coefs.col(1).array() += 0.5; + patches.addPatch(basis.makeGeometry(coefs)); + patches.computeTopology(); + // Create bases + gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) + + // Set degree + bases.setDegree( bases.maxCwiseDegree() + numElevate); + + // h-refine each basis + for (int r =0; r < numRefine; ++r) + bases.uniformRefine(); + numRefine = 0; + + gsInfo << "Patches: "<< patches.nPatches() <<", degree: "<< bases.minCwiseDegree() <<"\n"; + + // get from XML + real_t rho = 3000; + real_t E = 4000000; + real_t nu = 0.3; + + // get from XML ?? + // Set the interface for the precice coupling + std::vector > patchCouplingInterfaces(patches.nPatches()); + patchCouplingInterfaces[0].push_back(patchSide(0,boundary::east)); + patchCouplingInterfaces[0].push_back(patchSide(0,boundary::west)); + patchCouplingInterfaces[1].push_back(patchSide(1,boundary::east)); + patchCouplingInterfaces[1].push_back(patchSide(1,boundary::north)); + patchCouplingInterfaces[1].push_back(patchSide(1,boundary::west)); + + /* + * Initialize the preCICE participant + * + * + */ + // get from XML + std::string participantName = "Solid"; + gsPreCICE participant(participantName, precice_config); + + /* + * Data initialization + * + * This participant manages the geometry. The follow meshes and data are made available: + * + * - Meshes: + * + KnotMesh This mesh contains the knots as mesh vertices + * + ControlPointMesh: This mesh contains the control points as mesh vertices + * + ForceMesh: This mesh contains the integration points as mesh vertices + * + * - Data: + * + ControlPointData: This data is defined on the ControlPointMesh and stores the displacement of the control points + * + ForceData: This data is defined on the ForceMesh and stores pressure/forces + */ + + std::string SolidMesh = "Solid-Mesh"; + std::string StressData = "Stress"; + std::string DisplacementData = "Displacement"; + + // Get the quadrature nodes on the coupling interface + gsOptionList quadOptions = gsAssembler<>::defaultOptions(); + + // Get the quadrature points + + // Containers for uv and xy coordinates of the quadrature points per patch + std::vector> patchQuad_uv(patches.nPatches()); + std::vector> patchQuad_xy(patches.nPatches()); + // Global containers for uv and xy coordinates of the quadrature points (to be sent to preCICE) + gsMatrix<> quad_uv; + gsMatrix<> quad_xy; + + // Offset per patch + std::vector patchOffsets(patches.nPatches()+1); + patchOffsets[0] = 0; // first patch starts at 0 + for (size_t p=0; p!=patches.nPatches(); ++p) + { + // Get the quadrature points for the current patch + patchQuad_uv[p] = gsQuadrature::getAllNodes(bases.basis(p),quadOptions,patchCouplingInterfaces[p]); + patchQuad_xy[p] = patches.patch(p).eval(patchQuad_uv[p]); + patchOffsets[p+1] = patchOffsets[p]+patchQuad_uv[p].cols(); // store the number of quadrature points per patch + } + + quad_uv.resize(patches.domainDim(), patchOffsets.back()); + quad_xy.resize(patches.geoDim(), patchOffsets.back()); + for (size_t p=0; p!=patches.nPatches(); ++p) + { + // Copy the quadrature points of the current patch to the global containers + quad_uv.block(0,patchOffsets[p],patches.domainDim(),patchQuad_uv[p].cols()) = patchQuad_uv[p]; + quad_xy.block(0,patchOffsets[p],patches.geoDim(),patchQuad_xy[p].cols()) = patchQuad_xy[p]; + } + gsVector quad_xyIDs; // needed for writing + participant.addMesh(SolidMesh,quad_xy,quad_xyIDs); + + // Check if initial data is required (call before initialize) + bool needsInitialData = participant.requiresInitialData(); + + // Define precice interface + real_t precice_dt = participant.initialize(); + + // Write initial data after initialize if required + if (needsInitialData) + { + gsInfo << "Writing initial displacement data...\n"; + gsMatrix<> initialDisplacement(patches.geoDim(), quad_xy.cols()); + initialDisplacement.setZero(); // Initial displacement is zero + participant.writeData(SolidMesh, DisplacementData, quad_xyIDs, initialDisplacement); + } + +// ---------------------------------------------------------------------------------------------- + + // Define boundary conditions + gsBoundaryConditions<> bcInfo; + // Dirichlet side + gsConstantFunction<> g_D(0,patches.geoDim()); + // Coupling side + // gsFunctionExpr<> g_C("1","0",patches.geoDim()); + + // get from XML ??? + // Initialize quad_stress matrix for storing stress data from PreCICE + gsMatrix<> quad_stress(2,quad_xy.cols()); + quad_stress.setZero(); + + // Create lookup function with multiple patches + gsLookupFunction g_L(patches.nPatches()); + + // Add lookup function data for each patch separately + for (size_t p = 0; p != patches.nPatches(); ++p) + { + // Get stress data for this patch + gsMatrix<> patch_stress = quad_stress.block(0, patchOffsets[p], 2, patchQuad_xy[p].cols()); + g_L.add(patchQuad_xy[p], patch_stress); + } + + // gsPreCICEFunction g_C(&participant,SolidMesh,StressData,patches,patches.geoDim(),false); + // Add all BCs + // Coupling interface + for (size_t p=0; p!=patchCouplingInterfaces.size(); ++p) + for (size_t i=0; i!=patchCouplingInterfaces[p].size(); ++i) + bcInfo.addCondition(patchCouplingInterfaces[p][i].patch, patchCouplingInterfaces[p][i].side(), condition_type::neumann, &g_L, -1, true); + + // Bottom side (prescribed temp) + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &g_D, 0); + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &g_D, 1); + // Assign geometry map + bcInfo.setGeoMap(patches); + +// ---------------------------------------------------------------------------------------------- + // source function, rhs + gsConstantFunction<> g(0.,0.,2); + + // source function, rhs + gsConstantFunction<> gravity(0.,0.,2); + + // creating mass assembler + gsMassAssembler massAssembler(patches,bases,bcInfo,gravity); + massAssembler.options().setReal("Density",rho); + massAssembler.assemble(); + + // creating stiffness assembler. + gsElasticityAssembler assembler(patches,bases,bcInfo,g); + assembler.options().setReal("YoungsModulus",E); + assembler.options().setReal("PoissonsRatio",nu); + assembler.options().setInt("MaterialLaw",material_law::saint_venant_kirchhoff); + assembler.assemble(); + + gsMatrix Minv; + gsSparseMatrix<> M = massAssembler.matrix(); + gsSparseMatrix<> K = assembler.matrix(); + gsSparseMatrix<> K_T; + + // Time step + real_t dt = precice_dt; + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + // RHS of the projection + gsMatrix<> solVector; + solVector.setZero(assembler.numDofs(),1); + + std::vector > fixedDofs = assembler.allFixedDofs(); + // Assemble the RHS + gsVector<> U_checkpoint, V_checkpoint, A_checkpoint, U, V, A; + + U_checkpoint = U = gsVector::Zero(assembler.numDofs(),1); + V_checkpoint = V = gsVector::Zero(assembler.numDofs(),1); + A_checkpoint = A = gsVector::Zero(assembler.numDofs(),1); + + // Define the solution collection for Paraview + gsParaviewCollection collection("./output/solution"); + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&fixedDofs](gsMatrix const &x, gsSparseMatrix & m) { + // to do: add time dependency of forcing + assembler.assemble(x, fixedDofs); + m = assembler.matrix(); + return true; + }; + + // Function for the Residual + gsStructuralAnalysisOps::TResidual_t Residual = [&assembler,&fixedDofs](gsMatrix const &x, real_t t, gsVector & result) + { + assembler.assemble(x,fixedDofs); + result = assembler.rhs(); + return true; + }; + + // Function for the Residual + gsStructuralAnalysisOps::TForce_t TForce = [&assembler](real_t, gsVector & result) + { + assembler.assemble(); + result = assembler.rhs(); + return true; + }; + + + gsSparseMatrix<> C = gsSparseMatrix<>(assembler.numDofs(),assembler.numDofs()); + gsStructuralAnalysisOps::Damping_t Damping = [&C](const gsVector &, gsSparseMatrix & m) { m = C; return true; }; + gsStructuralAnalysisOps::Mass_t Mass = [&M]( gsSparseMatrix & m) { m = M; return true; }; + gsStructuralAnalysisOps::Stiffness_t Stiffness = [&K]( gsSparseMatrix & m) { m = K; return true; }; + + // if (nonlinear) + // timeIntegrator = new gsDynamicNewmark(Mass,Damping,Jacobian,Residual); + // else + // timeIntegrator = new gsDynamicNewmark(Mass,Damping,Stiffness,TForce); + + + + std::unique_ptr> timeIntegrator; + if (method==1) + timeIntegrator = std::make_unique>(Mass,Damping,Jacobian,Residual); + else if (method==2) + timeIntegrator = std::make_unique>(Mass,Damping,Jacobian,Residual); + else if (method==3 && nonlinear==false) + timeIntegrator = std::make_unique>(Mass,Damping,Stiffness,TForce); + else if (method==3 && nonlinear==true) + timeIntegrator = std::make_unique>(Mass,Damping,Jacobian,Residual); + else if (method==4) + timeIntegrator = std::make_unique>(Mass,Damping,Jacobian,Residual); + else if (method==5) + { + timeIntegrator = std::make_unique>(Mass,Damping,Jacobian,Residual); + timeIntegrator->options().setReal("gamma",1.4); + } + else if (method==6) + timeIntegrator = std::make_unique>(Mass,Damping,Jacobian,Residual); + else + GISMO_ERROR("Method "<options().setReal("DT",dt); + timeIntegrator->options().setReal("TolU",1e-6); + timeIntegrator->options().setSwitch("Verbose",true); + + real_t time = 0; + + // Plot initial solution + if (plot) + { + gsMultiPatch<> solution; + assembler.constructSolution(solVector,fixedDofs,solution); + + gsField<> solField(patches,solution); + std::string fileName = "./output/solution" + util::to_string(timestep); + gsWriteParaview<>(solField, fileName, 500); + for (size_t p=0; p!=patches.nPatches(); ++p) + { + fileName = "solution" + util::to_string(timestep) + util::to_string(p); + collection.addPart(fileName,time,"Solution",p); + } + } + + gsMatrix<> refPoints(2,1); + refPoints.col(0)<<0.5,1; + gsVector refPatches(1); + refPatches[0] = 1; // patch 0 + + gsStructuralAnalysisOutput writer("pointData.csv",refPoints); + writer.init({"x","y"},{"time"}); // point1 - x, point1 - y, time + + gsMatrix<> pointDataMatrix; + gsMatrix<> otherDataMatrix(1,1); + + // Time integration loop + gsInfo << "Starting coupling time loop...\n"; + while (participant.isCouplingOngoing()) + { + if (participant.requiresWritingCheckpoint()) + { + U_checkpoint = U; + V_checkpoint = V; + A_checkpoint = A; + + gsInfo<<"Checkpoint written:\n"; + gsInfo<<"\t ||U|| = "< patch_stress = quad_stress.block(0, patchOffsets[p], 2, patchQuad_xy[p].cols()); + g_L.set(p, patchQuad_xy[p], patch_stress); + } + + // solve gismo timestep + gsInfo << "Solving timestep " << time << "...\n"; + timeIntegrator->step(time,dt,U,V,A); + solVector = U; + gsInfo<<"Finished\n"; + + // potentially adjust non-matching timestep sizes + dt = std::min(dt,precice_dt); + + + gsMultiPatch<> solution; + assembler.constructSolution(solVector,fixedDofs,solution); + gsMatrix<> result(patches.geoDim(),patchOffsets.back()); + gsMatrix<> tmpResult; + for (size_t p=0; p!=patches.nPatches(); ++p) + { + // evaluate the solution at the quadrature points + solution.patch(p).eval_into(patchQuad_uv[p],tmpResult); + result.block(0,patchOffsets[p],patches.geoDim(),patchQuad_uv[p].cols()) = tmpResult; + } + participant.writeData(SolidMesh,DisplacementData,quad_xyIDs,result); + + // do the coupling + precice_dt =participant.advance(dt); + + if (participant.requiresReadingCheckpoint()) + { + U = U_checkpoint; + V = V_checkpoint; + A = A_checkpoint; + timestep = timestep_checkpoint; + } + else + { + // gsTimeIntegrator advances + // advance variables + time += dt; + timestep++; + + gsField<> solField(patches,solution); + if (timestep % plotmod==0) // Generate Paraview output for visualization + { + if (plot) + { + // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + std::string fileName = "./output/solution" + util::to_string(timestep); + gsWriteParaview<>(solField, fileName, 500); + for (size_t p=0; p!=patches.nPatches(); ++p) + { + fileName = "solution" + util::to_string(timestep) + util::to_string(p); + collection.addPart(fileName,time,"Solution",p); + } + } + + gsMatrix<> tmpData; + pointDataMatrix.resize(solution.targetDim(),refPoints.cols()); + for (size_t k=0; k!=refPoints.cols(); ++k) + { + // evaluate the solution at the reference points + solution.patch(refPatches[k]).eval_into(refPoints.col(k),tmpData); + if (pointDataMatrix.rows()==0) + pointDataMatrix.resize(tmpData.rows(),refPoints.cols()); + pointDataMatrix.col(k) = tmpData.col(0); + } + + solution.patch(0).eval_into(refPoints,pointDataMatrix); + otherDataMatrix< -#include -#include +#include + namespace gismo { +template +class gsLookupFunctionSingle; +/** + * @brief Multipatch lookup function evaluator + * @details This class enables lookup functions on multipatch geometries where each patch + * can have its own lookup function mapping points to data values + * @param T Number format + */ +template +class gsLookupFunction : public gsFunctionSet +{ +public: + typedef memory::shared_ptr< gsLookupFunction > Ptr; + typedef memory::unique_ptr< gsLookupFunction > uPtr; + typedef gsLookupFunctionSingle Piece_t; + + typedef std::vector Container; + + // Constructor + gsLookupFunction(size_t nPatches = 0) + : m_container() + { + m_container.reserve(nPatches); + } + /// Constructor with points and data for a single patch + gsLookupFunction(const gsMatrix & points, + const gsMatrix & data) + : gsLookupFunction(1) + { + this->add(points, data); + } + + /// Destructor + ~gsLookupFunction() + { + } + + GISMO_CLONE_FUNCTION(gsLookupFunction) + + /// Number of pieces + index_t nPieces() const override + { return m_container.size(); } + + /// Access a piece + const gsFunction & piece(const index_t k) const override + { + GISMO_ASSERT(k < m_container.size(), "Piece index out of bounds"); + return m_container[k]; + } + + /// See \a gsFunctionSet + index_t size() const override + { return m_container.size(); } + + /// Domain dimension + short_t domainDim() const override + { + if (m_container.size() > 0) + return m_container[0].domainDim(); + return 0; + } + + /// Target dimension + short_t targetDim() const override + { + if (m_container.size() > 0) + return m_container[0].targetDim(); + return 0; + } + + /// Evaluate function - requires patch index (no default evaluation across all patches) + void eval_into(const gsMatrix& u, gsMatrix& result) const override + { + GISMO_ERROR("gsLookupFunction requires patch index for evaluation. Use eval_into(patch, u, result)"); + } + + /// Evaluate derivatives - requires patch index (no default evaluation across all patches) + void deriv_into(const gsMatrix& u, gsMatrix& result) const override + { + GISMO_ERROR("gsLookupFunction requires patch index for evaluation. Use deriv_into(patch, u, result)"); + } + + /// Evaluate second derivatives - requires patch index (no default evaluation across all patches) + void deriv2_into(const gsMatrix& u, gsMatrix& result) const override + { + GISMO_ERROR("gsLookupFunction requires patch index for evaluation. Use deriv2_into(patch, u, result)"); + } + + /// Evaluate the function at points u for specific patch + void eval_into(const index_t patch, const gsMatrix& u, gsMatrix& result) const + { + GISMO_ASSERT(patch < m_container.size(), "Patch index out of bounds"); + m_container[patch].eval_into(u, result); + } + + /// Evaluate derivatives for specific patch + void deriv_into(const index_t patch, const gsMatrix& u, gsMatrix& result) const + { + GISMO_ASSERT(patch < m_container.size(), "Patch index out of bounds"); + m_container[patch].deriv_into(u, result); + } + + /// Evaluate second derivatives for specific patch + void deriv2_into(const index_t patch, const gsMatrix& u, gsMatrix& result) const + { + GISMO_ASSERT(patch < m_container.size(), "Patch index out of bounds"); + m_container[patch].deriv2_into(u, result); + } + + /// Get the lookup function container + const Container & container() const { return m_container; } + + /// Get the lookup function container (non-const) + Container & container() { return m_container; } + + /// Add lookup function (appends to the end) + void add(const gsMatrix & points, const gsMatrix & data) + { + m_container.push_back(gsLookupFunctionSingle(points, data)); + } + + /// Set lookup function for a specific patch index + void set(index_t patch, const gsMatrix & points, const gsMatrix & data) + { + GISMO_ASSERT(patch < m_container.size(), "Patch index out of bounds"); + m_container[patch] = gsLookupFunctionSingle(points, data); + } + +protected: + + Container m_container; +}; /** * @brief Class defining a function that looks up registered data on points. @@ -31,7 +162,7 @@ namespace gismo */ template -class gsLookupFunction : public gsFunction +class gsLookupFunctionSingle : public gsFunction { struct Compare @@ -46,16 +177,18 @@ class gsLookupFunction : public gsFunction typedef gsGeometry Base; /// Shared pointer for gsLookupFunction - typedef memory::shared_ptr< gsLookupFunction > Ptr; + typedef memory::shared_ptr< gsLookupFunctionSingle > Ptr; /// Unique pointer for gsLookupFunction - typedef memory::unique_ptr< gsLookupFunction > uPtr; + typedef memory::unique_ptr< gsLookupFunctionSingle > uPtr; - // /// Default constructor - // gsLookupFunction() { } + /// Default constructor + gsLookupFunctionSingle() { + // Note: m_points and m_data are empty, update() will be called when data is added + } /** - * @brief Constructs a new instance of the gsLookupFunction + * @brief Constructs a new instance of the gsLookupFunctionSingle * * @param interface The precice::SolverInterface (see \a gsPreCICE) * @param[in] meshName The ID of the mesh on which the data is located @@ -63,8 +196,8 @@ class gsLookupFunction : public gsFunction * @param[in] patches The geometry * @param[in] parametric Specifies whether the data is defined on the parametric domain or not */ - gsLookupFunction( const gsMatrix & points, - const gsMatrix & data ) + gsLookupFunctionSingle( const gsMatrix & points, + const gsMatrix & data ) : m_points(points), m_data(data) @@ -75,12 +208,12 @@ class gsLookupFunction : public gsFunction /// Constructs a function pointer static uPtr make( const gsMatrix & points, const gsMatrix & data ) - { return uPtr(new gsLookupFunction(points, data)); } + { return uPtr(new gsLookupFunctionSingle(points, data)); } - GISMO_CLONE_FUNCTION(gsLookupFunction) + GISMO_CLONE_FUNCTION(gsLookupFunctionSingle) /// Access a piece - const gsLookupFunction & piece(const index_t) const + const gsLookupFunctionSingle & piece(const index_t) const { return *this; } @@ -163,17 +296,18 @@ class gsLookupFunction : public gsFunction /// See \a gsFunction virtual std::ostream &print(std::ostream &os) const { - os << "gsLookupFunction\n"; + os << "gsLookupFunctionSingle\n"; return os; } protected: - const gsMatrix & m_points; - const gsMatrix & m_data; + gsMatrix m_points; + gsMatrix m_data; std::map,index_t,Compare> m_map; }; -} +} // namespace gismo +