diff --git a/examples/perpendicular-flap-shell-gismo.cpp b/examples/perpendicular-flap-shell-gismo.cpp new file mode 100644 index 0000000..5101d52 --- /dev/null +++ b/examples/perpendicular-flap-shell-gismo.cpp @@ -0,0 +1,460 @@ +/** + * @file perpendicular-flap-shell-gismo.cpp + * @brief Calculate mechanical response of perpendicular flap shell element + * + * This program calculates the mechanical response of a perpendicular flap shell element, + * including deformation, stress and strain. It uses the gsPreCICE library for structural + * analysis and gsThinShellAssembler for shell element analysis. + * + * For Kirchhoff–Love shells, a mapping is established between the mid-surface (P) + * and the top/bottom surfaces (Q). In this version, the mapping is updated every time step + * based on the current mid-surface geometry: The current unit normal is recalculated, + * and the corresponding thickness offset is computed for updated 3D control points. + * + * This ensures that even for large deformations and rotations the generated volume + * (i.e. the top/bottom surfaces) remain perpendicular to the current mid-surface. + * + * Author: J. Li, H. Verhelst + */ + +#include +#include +#include +#include +#include + +#include + +#ifdef gsKLShell_ENABLED +#include +#include +#endif + +#ifdef gsStructuralAnalysis_ENABLED +#include +#include +#include +#include +#include +#include +#include +#endif + +using namespace gismo; + +int main(int argc, char *argv[]) +{ + //////////////////////////////////////////////////////////////////////////// + // Parse command line options + bool plot = false; + bool write = false; + bool get_readTime = false; + bool get_writeTime = false; + bool nonlinear = false; + index_t plotmod = 1; + index_t numRefine = 1; + index_t numElevate = 1; + std::string precice_config; + index_t n = 5; + index_t m = 5; + index_t degree = 3; + int method = 3; // 1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson, 6: RK4 + std::string output(""); + std::string dirname = "./output"; + + 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", "method", "1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson", method); + cmd.addSwitch("write", "Create a file with point data", write); + cmd.addSwitch("readTime", "Get the read time", get_readTime); + cmd.addSwitch("writeTime", "Get the write time", get_writeTime); + cmd.addInt("n", "dof1", "Number of basis function in one direction", n); + cmd.addInt("d", "degree", "Degree of a surface", degree); + cmd.addString("o", "output", "Name of the output file.", output); + try { cmd.getValues(argc, argv); } catch (int rv) { return rv; } + + GISMO_ASSERT(gsFileManager::fileExists(precice_config), "No precice config file has been defined"); + + // Adjust values if necessary + degree = math::max((index_t)(0), degree); + n = math::max(n, degree + 1); + m = math::max(m, degree + 1); + + gsInfo << "----------------------\n\n" + << "n: " << n << "\n" + << "degree: " << degree << "\n" + << "output: " << output << "\n" + << "----------------------\n\n"; + + //////////////////////////////////////////////////////////////////////////// + // Construct the 2D mid-surface geometry + gsKnotVector<> kv1(0, 1, 0, 2); + gsKnotVector<> kv2(0, 1, n - degree - 1, degree + 1); + gsTensorBSplineBasis<2, real_t> basis(kv1, kv2); + gsMatrix<> greville = basis.anchors(); + + // Generate control points for the vertical flap. + // Flap is positioned in the YZ plane with x = 0. + gsMatrix<> coefs(greville.cols(), 3); + for (index_t col = 0; col != greville.cols(); col++) + { + real_t x = greville(0, col); + real_t y = greville(1, col); + coefs(col, 0) = 0; + coefs(col, 1) = y; + coefs(col, 2) = x; + } + gsTensorBSpline<2, real_t> surface(basis, coefs); + gsExprEvaluator<> ev; + ev.setIntegrationElements(basis); + auto G = ev.getMap(surface); + + //////////////////////////////////////////////////////////////////////////// + // Define thickness for shell and create 3D volume for coupling interface + gsFunctionExpr<> thickness("0.1", 2); + real_t shell_thickness = 0.1; + + // Create 3D volume using surfaceToVolume for coupling interface only + gsTensorBSpline<3, real_t> volume3D = surfaceToVolume(surface, thickness); + gsGeometry<>::uPtr volume = memory::make_unique(new gsTensorBSpline<3, real_t>(volume3D)); + const gsTensorBSplineBasis<3, real_t>& vbasis = volume3D.basis(); + + //////////////////////////////////////////////////////////////////////////// + // Use the original 2D surface as the mid-surface geometry for shell analysis + gsMultiPatch<> mid_surface_geom; + mid_surface_geom.addPatch(surface.clone()); + gsWriteParaview(mid_surface_geom, "surf_mid", 1000, true); + + gsOptionList quadOptions = gsExprAssembler<>::defaultOptions(); + gsMatrix<> quad_shell_uv = gsQuadrature::getAllNodes(basis.basis(0), quadOptions); + gsMatrix<> quad_shell_xy = mid_surface_geom.patch(0).eval(quad_shell_uv); + + // Embed 2D geometry into 3D for coupling. + gsMultiPatch<> solutions; + gsConstantFunction<> g(0., 0., 3); + gsConstantFunction<> gravity(0., 0., 3); + + //////////////////////////////////////////////////////////////////////////// + // PreCICE coupling initialization. + std::string participantName = "Solid"; + gsPreCICE participant(participantName, precice_config); + + + std::string SolidMesh = "Solid-Mesh"; + std::string StressData = "Stress"; + std::string DisplacementData = "Displacement"; + + // For coupling, we get front and back faces separately to ensure correct ordering + std::vector frontInterface, backInterface; + frontInterface.push_back(patchSide(0, boundary::front)); // Top surface (w=1) + backInterface.push_back(patchSide(0, boundary::back)); // Bottom surface (w=0) + + + // Get 3D quadrature points separately for each surface + gsMatrix<> quad_uv_front = gsQuadrature::getAllNodes(vbasis.basis(0), quadOptions, frontInterface); + gsMatrix<> quad_uv_back = gsQuadrature::getAllNodes(vbasis.basis(0), quadOptions, backInterface); + + + // Evaluate positions + gsMatrix<> quad_xy_front = volume->eval(quad_uv_front); + gsMatrix<> quad_xy_back = volume->eval(quad_uv_back); + + // Combine them in a known order: first all front points, then all back points + gsMatrix<> quad_xy(3, quad_xy_front.cols() + quad_xy_back.cols()); + quad_xy.leftCols(quad_xy_front.cols()) = quad_xy_front; + quad_xy.rightCols(quad_xy_back.cols()) = quad_xy_back; + + // Also combine parametric coordinates for later use + gsMatrix<> quad_uv(3, quad_uv_front.cols() + quad_uv_back.cols()); + quad_uv.leftCols(quad_uv_front.cols()) = quad_uv_front; + quad_uv.rightCols(quad_uv_back.cols()) = quad_uv_back; + + gsWriteParaviewPoints(quad_xy, "quadPointsAll"); + + gsVector quadPointIDs; + participant.addMesh(SolidMesh, quad_xy, quadPointIDs); + + real_t precice_dt = participant.initialize(); + + gsBoundaryConditions<> bcInfo; + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, 0, 0, false, 0); + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, 0, 0, false, 1); + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, 0, 0, false, 2); + bcInfo.addCondition(0, boundary::south, condition_type::clamped, 0, 0, false, 0); + bcInfo.addCondition(0, boundary::east, condition_type::clamped, 0, 0, false, 2); + bcInfo.addCondition(0, boundary::west, condition_type::clamped, 0, 0, false, 2); + bcInfo.setGeoMap(mid_surface_geom); + + //////////////////////////////////////////////////////////////////////////// + // Material properties and assembler setup. + real_t rho = 3000; + real_t E = 4e6; + real_t nu = 0.3; + gsFunctionExpr<> E_modulus(std::to_string(E), 3); + gsFunctionExpr<> PoissonRatio(std::to_string(nu), 3); + gsFunctionExpr<> Density(std::to_string(rho), 3); + gsFunctionExpr<> thickness_3d("0.1", 3); + + std::vector*> parameters(2); + parameters[0] = &E_modulus; + parameters[1] = &PoissonRatio; + + gsOptionList options; + // Initialize force data as zero for all shell quadrature points + // Forces will only be non-zero at specific boundary points + gsMatrix<> quadPointsData(3, quad_shell_xy.cols()); + quadPointsData.setZero(); + + // Create lookup function for shell forces + gsLookupFunction surfForce(quad_shell_xy, quadPointsData); + + gsMaterialMatrixBase::uPtr materialMatrix; + options.addInt("Material", "Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden", 0); + options.addInt("Implementation", "Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral", 1); + materialMatrix = getMaterialMatrix<3, real_t>(mid_surface_geom, thickness_3d, parameters, Density, options); + + gsMultiBasis<> bases(basis); + gsThinShellAssemblerBase* assembler; + // gsFunctionExpr dummysurf("t","0","0", 3); + // dummysurf.set_t(0); // Set time to zero for initial conditions + assembler = new gsThinShellAssembler<3, real_t, true>(mid_surface_geom, bases, bcInfo, surfForce, materialMatrix); // Use linear for now + gsOptionList assemblerOptions = options.wrapIntoGroup("Assembler"); + assembler->assemble(); + assembler->setOptions(assemblerOptions); + + // Assemble constant mass and stiffness matrices. + assembler->assembleMass(); + gsSparseMatrix<> M = assembler->massMatrix(); + gsInfo << "Mass norm" << M.norm() << "\n"; + assembler->assemble(); + gsSparseMatrix<> K = assembler->matrix(); + + gsFileManager::mkdir(dirname); + gsParaviewCollection collection(dirname + "/solution"); + gsParaviewCollection collection_surface(dirname + "/deformed_surface"); + gsParaviewCollection collection_volume(dirname + "/deformed_volume"); + + //////////////////////////////////////////////////////////////////////////// + // Time stepping and dynamic solver setup. + real_t t_read = 0, t_write = 0; + real_t dt = precice_dt; + + gsStructuralAnalysisOps::Jacobian_t Jacobian = + [&assembler, &solutions](gsMatrix const &x, gsSparseMatrix &m) + { + assembler->constructSolution(x, solutions); + ThinShellAssemblerStatus status = assembler->assembleMatrix(solutions); + m = assembler->matrix(); + return true; + }; + + gsStructuralAnalysisOps::TResidual_t Residual = + [&assembler, &solutions](gsMatrix const &x, real_t t, gsVector &result) + { + assembler->constructSolution(x, solutions); + ThinShellAssemblerStatus status = assembler->assembleVector(solutions); + 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(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; }; + + + gsDynamicBase* timeIntegrator; + if (method == 1) + timeIntegrator = new gsDynamicExplicitEuler(Mass, Damping, Jacobian, Residual); + else if (method == 2) + timeIntegrator = new gsDynamicImplicitEuler(Mass, Damping, Jacobian, Residual); + else if (method==3 && nonlinear==false) + timeIntegrator = new gsDynamicNewmark(Mass,Damping,Stiffness,TForce); + else if (method==3 && nonlinear==true) + timeIntegrator = new gsDynamicNewmark(Mass,Damping,Jacobian,Residual); + else if (method == 4) + timeIntegrator = new gsDynamicBathe(Mass, Damping, Jacobian, Residual); + else if (method == 5) + { + timeIntegrator = new gsDynamicWilson(Mass, Damping, Jacobian, Residual); + timeIntegrator->options().setReal("gamma", 1.4); + } + else if (method == 6) + timeIntegrator = new gsDynamicRK4(Mass, Damping, Jacobian, Residual); + else + GISMO_ERROR("Method " << method << " not known"); + + timeIntegrator->options().setReal("DT", dt); + timeIntegrator->options().setReal("TolU", 1e-3); + timeIntegrator->options().setSwitch("Verbose", true); + + // Initialize displacement, velocity and acceleration (all zero). + gsMatrix<> solVector = gsMatrix::Zero(assembler->numDofs(), 1); + gsVector<> F = assembler->rhs(); + + gsVector<> F_checkpoint, U_checkpoint, V_checkpoint, A_checkpoint, U, V, A; + F_checkpoint = F; + U_checkpoint = U = gsVector::Zero(assembler->numDofs(), 1); + V_checkpoint = V = gsVector::Zero(assembler->numDofs(), 1); + A_checkpoint = A = gsVector::Zero(assembler->numDofs(), 1); + + index_t timestep = 0, timestep_checkpoint = 0; + real_t time = 0; + + //////////////////////////////////////////////////////////////////////////// + // Main time stepping loop. + while (participant.isCouplingOngoing()) + { + if (participant.requiresWritingCheckpoint()) + { + U_checkpoint = U; + V_checkpoint = V; + A_checkpoint = A; + gsInfo << "Checkpoint written:\n" + << "\t ||U|| = " << U.norm() << "\n" + << "\t ||V|| = " << V.norm() << "\n" + << "\t ||A|| = " << A.norm() << "\n"; + timestep_checkpoint = timestep; + } + + // Read stress data from preCICE (3D forces on top and bottom surfaces) + gsMatrix<> stressDataRead; + participant.readData(SolidMesh, StressData, quadPointIDs, stressDataRead); + + // Now we know the exact ordering: first quad_xy_front.cols() points are from front, + // next quad_xy_back.cols() points are from back + index_t numFrontPoints = quad_xy_front.cols(); + index_t numBackPoints = quad_xy_back.cols(); + + + // Average forces from front (top) and back (bottom) surfaces to get mid-surface forces + gsMatrix<> avgForces(3, numFrontPoints); + + for (index_t i = 0; i < numFrontPoints; ++i) + { + avgForces.col(i) = (stressDataRead.col(i) + stressDataRead.col(numFrontPoints + i)); + } + + // Extract 2D parametric coordinates from 3D front boundary + gsMatrix<> surf_quad_uv(2, numFrontPoints); + for (index_t i = 0; i < numFrontPoints; ++i) + { + surf_quad_uv(0, i) = quad_uv_front(0, i); // u parameter + surf_quad_uv(1, i) = quad_uv_front(1, i); // v parameter + } + + // Evaluate 2D positions for these parametric points + gsMatrix<> surf_quad_xy = surface.eval(surf_quad_uv); + // Update the quadPointsData with the average forces for the look up function + quadPointsData = avgForces; + // quadPointsData.setOnes(); + + if (get_readTime) + t_read += participant.readTime(); + assembler->assemble(); + F = assembler->rhs(); + + gsInfo << "Solving timestep " << time << "...\n"; + timeIntegrator->step(time, dt, U, V, A); + solVector = U; + gsInfo << "Finished\n"; + + dt = std::min(dt, precice_dt); + + // Construct the deformed solution of the mid-surface from current displacement U (updated mid-surface geometry) + gsMultiPatch<> solution = assembler->constructSolution(U); + // Calculate mid-surface displacement at integration points and control points (for output only, not used for new normal calculation) + gsMatrix<> MidPointDisp = solution.patch(0).eval(quad_shell_uv); + + // Calculate deformed surface and create deformed volume + gsTensorBSpline<2, real_t> deformed_surface(basis, solution.patch(0).coefs()); + + // Create deformed 3D volume using surfaceToVolume + gsTensorBSpline<3, real_t> deformed_volume3D = surfaceToVolume(deformed_surface, thickness); + + // Evaluate deformed positions at coupling points + gsMatrix<> deformed_quad_xy = deformed_volume3D.eval(quad_uv); + + // Calculate 3D displacement + gsMatrix<> displacementData = deformed_quad_xy - quad_xy; + + // Write 3D displacement data for coupling + participant.writeData(SolidMesh, DisplacementData, quadPointIDs, displacementData); + + if (get_writeTime) + t_write += participant.writeTime(); + + // Perform PreCICE coupling + precice_dt = participant.advance(dt); + + if (participant.requiresReadingCheckpoint()) + { + U = U_checkpoint; + V = V_checkpoint; + A = A_checkpoint; + timestep = timestep_checkpoint; + } + else + { + time += dt; + timestep++; + + assembler->constructDisplacement(U, solution); + + gsField<> solField(mid_surface_geom, solution); + if (timestep % plotmod == 0 && plot) + { + std::string fileName = dirname + "/solution" + util::to_string(timestep); + gsWriteParaview<>(solField, fileName, 500); + std::string fileNameSol = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileNameSol, time, ".vts"); + + // For 2D case, save the deformed surface instead of volume + std::string fileName2 = dirname + "/deformed_surface" + util::to_string(timestep); + gsWriteParaview(deformed_surface, fileName2, 1000, true); + std::string fileNameSurf = "deformed_surface" + util::to_string(timestep); + collection_surface.addTimestep(fileNameSurf, time, ".vts"); + + std::string fileName3 = dirname + "/deformed_volume" + util::to_string(timestep); + gsWriteParaview(deformed_volume3D, fileName3, 1000, true); + std::string fileNameVol = "deformed_volume" + util::to_string(timestep); + collection_volume.addTimestep(fileNameVol, time, ".vts"); + } + + } + } // end while coupling loop + + if (get_readTime) + gsInfo << "Read time: " << t_read << "\n"; + if (get_writeTime) + gsInfo << "Write time: " << t_write << "\n"; + + if (plot) + { + collection.save(); + collection_surface.save(); + collection_volume.save(); + + } + + delete assembler; + delete timeIntegrator; + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/examples/perpendicular-flap-vertex-gismo-volume.cpp b/examples/perpendicular-flap-vertex-gismo-volume.cpp new file mode 100644 index 0000000..c48cf89 --- /dev/null +++ b/examples/perpendicular-flap-vertex-gismo-volume.cpp @@ -0,0 +1,373 @@ +/** @file solid-gismo-elasticity.cpp + + @brief Elasticity participant for the PreCICE example "perpendicular-flap". + + + 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 + +#ifdef gsElasticity_ENABLED +#include +#include +#endif + +#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,24,2) ; + gsKnotVector<> KV3 (0,1,0,2) ; + + gsTensorBSplineBasis<3> basis(KV1,KV2,KV3); + gsMatrix<> coefs = basis.anchors(); + coefs.transposeInPlace(); + coefs.col(0).array() *= 0.1; + coefs.col(0).array() -= 0.05; + + + patches.addPatch(basis.makeGeometry(give(coefs))); + + // 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 couplingInterfaces(3); + couplingInterfaces[0] = patchSide(0,boundary::east); + couplingInterfaces[1] = patchSide(0,boundary::north); + couplingInterfaces[2] = patchSide(0,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 + gsMatrix<> quad_uv = gsQuadrature::getAllNodes(bases.basis(0),quadOptions,couplingInterfaces); // Quadrature points in the parametric domain + gsMatrix<> quad_xy = patches.patch(0).eval(quad_uv); // Quadrature points in the physical domain + gsVector quad_xyIDs; // needed for writing + participant.addMesh(SolidMesh,quad_xy,quad_xyIDs); + + // Define precice interface + real_t precice_dt = participant.initialize(); + +// ---------------------------------------------------------------------------------------------- + + // 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 ??? + gsMatrix<> quad_stress(3,quad_xy.cols()); + gsLookupFunction g_L(quad_xy,quad_stress); + // gsPreCICEFunction g_C(&participant,SolidMesh,StressData,patches,patches.geoDim(),false); + // Add all BCs + // Coupling interface + bcInfo.addCondition(0, boundary::north, condition_type::neumann , &g_L,-1,true); + bcInfo.addCondition(0, boundary::east, condition_type::neumann , &g_L,-1,true); + bcInfo.addCondition(0, boundary::west, 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); + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &g_D, 2); + + bcInfo.addCondition(0, boundary::front, condition_type::dirichlet, &g_D, 2); + bcInfo.addCondition(0, boundary::back, condition_type::dirichlet, &g_D, 2); + // Assign geometry map + bcInfo.setGeoMap(patches); + +// ---------------------------------------------------------------------------------------------- + // source function, rhs + gsConstantFunction<> g(0.,0., 0., 3); + + // source function, rhs + gsConstantFunction<> gravity(0.,0.,0.,3); + // 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; }; + + gsDynamicBase * timeIntegrator; + if (method==1) + timeIntegrator = new gsDynamicExplicitEuler(Mass,Damping,Jacobian,Residual); + else if (method==2) + timeIntegrator = new gsDynamicImplicitEuler(Mass,Damping,Jacobian,Residual); + else if (method==3 && nonlinear==false) + timeIntegrator = new gsDynamicNewmark(Mass,Damping,Stiffness,TForce); + else if (method==3 && nonlinear==true) + timeIntegrator = new gsDynamicNewmark(Mass,Damping,Jacobian,Residual); + else if (method==4) + timeIntegrator = new gsDynamicBathe(Mass,Damping,Jacobian,Residual); + else if (method==5) + { + timeIntegrator = new gsDynamicWilson(Mass,Damping,Jacobian,Residual); + timeIntegrator->options().setReal("gamma",1.4); + } + else if (method==6) + timeIntegrator = new gsDynamicRK4(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); + fileName = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileName,time,".vts"); + } + + gsMatrix<> points(3,1); + points.col(0)<<0.5,1,0; + + gsStructuralAnalysisOutput writer("pointData.csv",points); + writer.init({"x","y","z"},{"time"}); // point1 - x, point1 - y, point1 - z, time + + gsMatrix<> pointDataMatrix; + gsMatrix<> otherDataMatrix(1,1); + + // Time integration loop + while (participant.isCouplingOngoing()) + { + if (participant.requiresWritingCheckpoint()) + { + U_checkpoint = U; + V_checkpoint = V; + A_checkpoint = A; + + gsInfo<<"Checkpoint written:\n"; + gsInfo<<"\t ||U|| = "<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); + // write heat fluxes to interface + gsMatrix<> result(patches.geoDim(),quad_uv.cols()); + solution.patch(0).eval_into(quad_uv,result); + 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); + fileName = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileName,time,".vts"); + } + + } + } + } + + if (plot) + { + collection.save(); + } + + delete timeIntegrator; + return EXIT_SUCCESS; +} diff --git a/examples/perpendicular-flap/precice-config-shell-3d.xml b/examples/perpendicular-flap/precice-config-shell-3d.xml new file mode 100644 index 0000000..9a95f80 --- /dev/null +++ b/examples/perpendicular-flap/precice-config-shell-3d.xml @@ -0,0 +1,64 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/gsPreCICEUtils.h b/gsPreCICEUtils.h index 334b659..4a7ff60 100644 --- a/gsPreCICEUtils.h +++ b/gsPreCICEUtils.h @@ -385,4 +385,124 @@ typename gsBasis::Ptr knotMatrixToBasis(const gsMatrix & knots) return memory::make_shared_not_owned(basis); } +/** + * @brief Creates a 3D volume from a 2D surface by extruding along its normals with variable thickness. + * + * This function takes a 2D tensor B-spline surface and extrudes it according to a scalar thickness function + * to create a 3D tensor B-spline volume. The extrusion is performed along the surface normals, with half + * of the thickness extending in each direction (positive and negative). + * + * @tparam T Numeric type (e.g., real_t, double) + * @param surface Input tensor B-spline surface (2D parametric domain) + * @param thickness Function that maps from the surface parameter space to thickness values + * + * @return A 3D tensor B-spline volume created by extruding the input surface + * + * @note The function supports both 2D and 3D target dimensions for the input surface. + * For 2D surfaces, the normal vector is assumed to be (0,0,1). + * For 3D surfaces, proper surface normals are computed. + * + * @throws AssertionError if the thickness function is not defined on a 2D domain + * @throws AssertionError if the thickness function is not scalar-valued + * @throws AssertionError if the surface is not defined on a 2D domain + * @throws Error if the surface's target dimension is not 2D or 3D + */ +template +gsTensorBSpline<3, T> surfaceToVolume(const gsTensorBSpline<2, T> & surface, const gsFunction & thickness) +{ + // Validate input parameters + GISMO_ASSERT(thickness.domainDim() == 2,"Thickness function must be defined on a 2D domain."); + GISMO_ASSERT(thickness.targetDim() == 1,"Thickness function must be scalar valued."); + GISMO_ASSERT(surface.domainDim() == 2,"Surface must be defined on a 2D domain."); + + // Create a knot vector for the third dimension (thickness direction) + // Linear basis with 2 elements: [0,0,0.5,1,1] (degree 1, 2 internal knots) + gsKnotVector kvZ(0, 1, 1, 2); + + // Create a 3D tensor B-spline basis by combining the existing 2D basis with the new Z knot vector + gsTensorBSplineBasis<3, T> basis(surface.basis().knots(0), + surface.basis().knots(1), + kvZ); + + // Calculate control point count for the new volume + const index_t nCoefs = surface.coefs().rows(); + + // Allocate control points matrix for the volume (3 layers of control points) + gsMatrix coefs(3*nCoefs, surface.coefs().cols()); + + // Get parametric coordinates (anchors) of the surface control points + gsMatrix anchors = surface.basis().anchors(); + + // Matrix to store thickness values at each anchor point + gsMatrix thicknessValues; + + // Evaluate the thickness function at each anchor point + thickness.eval_into(anchors, thicknessValues); + + // Matrix to store surface normals + gsMatrix normals; + + // Handle different dimension cases for the surface + if (surface.targetDim() == 2) + { + // For 2D surfaces, use (0,0,1) as the normal vector + normals = gsMatrix::Zero(3, anchors.cols()); + normals.row(2).setConstant(1.0); // Set the z-component of the normal to 1 + } + else if (surface.targetDim() == 3) + { + // For 3D surfaces, compute proper surface normals at each anchor point + gsMapData md(NEED_NORMAL); + md.points = anchors; + surface.computeMap(md); + std::swap(md.normals, normals); // Get the computed normal vectors + } + else + { + GISMO_ERROR("The surface must be 2D or 3D."); + } + + // For each control point of the surface + for (index_t k = 0; k != anchors.cols(); k++) + { + // Normalize the normal vector + const gsVector & normal = normals.col(k).normalized(); + + // Get thickness value at this control point + const T & t = thicknessValues(0,k); + + // Create three layers of control points: + // Bottom layer: offset by -0.5*thickness along the normal + coefs.row(k) = surface.coefs().row(k) - 0.5*t*normal.transpose(); + // Middle layer: original surface control point + coefs.row(k + nCoefs) = surface.coefs().row(k); + // Top layer: offset by +0.5*thickness along the normal + coefs.row(k + 2*nCoefs) = surface.coefs().row(k) + 0.5*t*normal.transpose(); + } + + // Create and return the 3D tensor B-spline volume + return gsTensorBSpline<3, T>(basis, coefs); +} + +/** + * @brief Converts a 2D tensor B-spline surface to a 3D tensor B-spline volume. + * + * This function creates a volumetric representation of a surface by extruding the + * surface in the normal direction with a constant thickness. + * + * @tparam T Type for coordinates (typically float or double) + * + * @param surface The input 2D tensor B-spline surface + * @param thickness The constant thickness value for the extrusion (default = 1.0) + * + * @return A 3D tensor B-spline volume representing the extruded surface + * + * @see toVolume(const gsTensorBSpline<2, T> &, const gsFunction &) + */ +template +gsTensorBSpline<3, T> toVolume(const gsTensorBSpline<2, T> & surface, const T & thickness = 1.0) +{ + return surfaceToVolume(surface, gsConstantFunction(thickness, surface.basis().domainDim())); +} + } //namespace gismo