diff --git a/examples/partitioned-heat-conduction-IGA-IGA/.DS_Store b/examples/partitioned-heat-conduction-IGA-IGA/.DS_Store new file mode 100644 index 0000000..adfa1e4 Binary files /dev/null and b/examples/partitioned-heat-conduction-IGA-IGA/.DS_Store differ diff --git a/examples/partitioned-heat-conduction-IGA-IGA/dirichlet/partitioned-heat-conduction-IGA-Dirichlet.cpp b/examples/partitioned-heat-conduction-IGA-IGA/dirichlet/partitioned-heat-conduction-IGA-Dirichlet.cpp new file mode 100644 index 0000000..3e5f41b --- /dev/null +++ b/examples/partitioned-heat-conduction-IGA-IGA/dirichlet/partitioned-heat-conduction-IGA-Dirichlet.cpp @@ -0,0 +1,391 @@ +/** @file heat-equation-coupling.cpp + + @brief Heat equation participant for a double coupled heat equation + + 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): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) +*/ + + + +#include +#include +#include +#include + +using namespace gismo; + +int main(int argc, char *argv[]) +{ + //! [Parse command line] + bool plot = true; + index_t plotmod = 1; + index_t numRefine = 2; + index_t numElevate = 0; + short_t side = 0; + real_t alpha = 3; + real_t beta = 1.2; + real_t time = 0; + real_t k_temp = 1; + + std::string precice_config("../precice_config.xml"); + + gsCmdLine cmd("Flow over heated plate for PreCICE."); + cmd.addString( "c", "config", "PreCICE config file", precice_config ); + cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot); + // cmd.addInt("l","loadCase", "Load case: 0=constant load, 1='spring' load", loadCase); + cmd.addInt("m","plotmod", "Modulo for plotting, i.e. if plotmod==1, plots every timestep", plotmod); + //cmd.addSwitch("readTime", "Get the read time", get_readTime); + //cmd.addSwitch("writeTime", "Get the write time", get_writeTime); + 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"); + + /* + * Initialize the preCICE participant + * + * + * Dirichlet side: + * Receive: temperature change as Dirichlet boundary condition + * Write: flux vector q(u) to the other side as Neumann condition + * + */ + std::string participantName = "Dirichlet"; + gsPreCICE participant(participantName, precice_config); + + + + //! [Read input file] + + /* + * Data initialization + * + * This participant receives mesh information (knot vector, control points) from the Neumann, + * and it creates its own mesh based on the information. + * And writes heat flux, reads temperature from the Neumann participant. + * 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 + * + FluxMesh: This mesh contains the integration points as mesh vertices + * + * - Data: + * + ControlPointData: This data is defined on the ControlPointMesh and stores the temperature of the control points + * + FluzData: This data is defined on the ForceMesh and stores pressure/forces + */ + + gsMultiPatch<> patches; + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,0.0,1.0,1.0)); + for (int r =0; r < numRefine; ++r) + patches.uniformRefine(); + + gsMultiBasis<> basesDirichlet(patches); + // ---------------------------------------------------------------------------------------------- + + + + /// Create from patches and boundary/interface information + + //Mesh provided by Neumann + std::string GeometryKnotMesh = "Geometry-Knot-Mesh"; + std::string GeometryControlPointMesh = "Geometry-Control-Point-Mesh"; + + // Mesh provided by Dirichlet + std::string FluxKnotMesh = "Flux-Knot-Mesh"; + std::string FluxControlPointMesh = "Flux-Control-Point-Mesh"; + + // Data provided by Neumann + std::string TemperatureData = "Temperature-Data"; + + // Data provided by Dirichlet + std::string FluxControlPointData = "Flux-Control-Point-Data"; + + // Setup bounding box onto the force mesh + gsMatrix<> bbox(2,2); + bbox << -1e300, 1e300, // X dimension limits + -1e300, 1e300; // Y dimension limits + bbox.transposeInPlace(); + bbox.resize(1,bbox.rows()*bbox.cols()); + participant.setMeshAccessRegion(GeometryControlPointMesh, bbox); + + + // ---------------------------------------------------------------------------------------------- + + std::vector couplingInterface(1); + couplingInterface[0] = patchSide(0,boundary::north); + std::vector::uPtr> boundaries(couplingInterface.size()); + + + + + gsMultiBasis<> fluxBases; + std::vector::uPtr> couplingBoundaries(couplingInterface.size()); + + couplingBoundaries[0] = patches.patch(0).boundary(couplingInterface[0].side()); + + // for(index_t i= 0; i < couplingInterface.size();++i) + // { + // couplingBoundaries[i] = patches.patch(0).boundary(couplingInterface[i].side()); // Add boundary coefficients to a vector + // auto& basis_temp = couplingBoundaries[i]->basis(); // Assuming this returns a pointer to gsBasis + // fluxBases.addBasis(&basis_temp); + // } + + // Add flux knot mesh + gsVector<> fluxKnotMatrix = knotsToVector(couplingBoundaries[0]->basis()); + gsDebugVar(fluxKnotMatrix); + participant.addMesh(FluxKnotMesh, fluxKnotMatrix.transpose()); + + + // Add flux control points mesh + gsVector fluxControlPointsIDs; + gsMatrix<> fluxControlPoints = couplingBoundaries[0]->coefs(); + gsDebugVar(fluxControlPoints); + participant.addMesh(FluxControlPointMesh,fluxControlPoints.transpose(), fluxControlPointsIDs); + + + real_t precice_dt = participant.initialize(); + + //Get the temperature mesh from direct-access="true"direct-access="true"the API + gsVector geometryKnotIDs; + gsMatrix<> geometryKnots; + participant.getMeshVertexIDsAndCoordinates(GeometryKnotMesh,geometryKnotIDs,geometryKnots); + + gsDebugVar(geometryKnots); + + //Get the temperature mesh from direct-access="true"direct-access="true"the API + gsVector geometryControlPointIDs; + gsMatrix<> geometryControlPoint; + participant.getMeshVertexIDsAndCoordinates(GeometryControlPointMesh,geometryControlPointIDs,geometryControlPoint); + + gsMatrix<> tempData(geometryControlPoint.rows(),geometryControlPoint.cols()); + tempData.setZero(); + + gsMultiPatch<> tempMesh; + gsBasis<> * basis = knotMatrixToBasis(geometryKnots.row(0)).get(); + tempMesh.addPatch(give(basis->makeGeometry(tempData.row(0).transpose()))); + gsDebugVar(tempMesh); + + + // Set external heat-flux to zero + gsConstantFunction<> f(beta-2-2*alpha,2); + gsFunctionExpr<> u_ex("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); + + gsBoundaryConditions<> bcInfo; + + // gsPreCICEFunction g_CD(&interface,meshName,(side==0 ? tempName : fluxName),patches,1); + gsFunction<> * g_C = &u_ex; + bcInfo.addCondition(0, boundary::north, condition_type::dirichlet, g_C, 0, false, 0); //For initialization, will be changed later + bcInfo.addCondition(0, boundary::west, condition_type::dirichlet , &u_ex, 0, false, 0); + bcInfo.addCondition(0, boundary::east, condition_type::dirichlet , &u_ex, 0, false, 0); + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet , &u_ex, 0, false, 0); + bcInfo.setGeoMap(patches); + + // ---------------------------------------------------------------------------------------------- + typedef gsExprAssembler<>::geometryMap geometryMap; + typedef gsExprAssembler<>::space space; + typedef gsExprAssembler<>::solution solution; + //-----------------------Get discretization space--------------- + gsExprAssembler<> A(1,1); + + gsInfo<<"Active options:\n"<< A.options() <<"\n"; + + A.setIntegrationElements(basesDirichlet); + + gsExprEvaluator<> ev(A); + + geometryMap G = A.getMap(patches); + + // Set the discretization space + space u = A.getSpace(basesDirichlet); // Use the Dirichlet basis function to discretize solution vector u + // u_h(x) = \sum_{i=1}^{N}u_i\varphi_i(x), where \varphi_i(x) is the basis function and u_i is the corresponding node values + + //----------------------RHS of the heat equation \frac{\partial u}{\partial t} = \nabla \cdot (\kappa \nabla u) + f --------------------- + // Set the source term + auto ff = A.getCoeff(f, G); + + // Set the solution and initial condition + gsMatrix<> solVector, solVector_ini; + solution u_sol = A.getSolution(u, solVector); + solution u_ini = A.getSolution(u, solVector); + + + // Initialize the system and assemble the mass matrix + // $M = \int_{\Omega} \varphi_i \varphi_j , d\Omega$ + u.setup(bcInfo, dirichlet::homogeneous, 0); + A.initSystem(); + gsDebugVar(A.numDofs()); + A.assemble( u * u.tr() * meas(G)); + gsSparseMatrix<> M = A.matrix(); + + // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner + gsSparseSolver<>::CGDiagonal solver; + + real_t dt = 0.1; + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + // $u_{\text{ini}}(x) = \int_{\Omega} \varphi_i u_{\text{ex}}(x) , d\Omega$ + auto uex = A.getCoeff(u_ex, G); + // RHS of the projection + u.setup(bcInfo, dirichlet::l2Projection, 0); + A.initSystem(); + A.assemble( u * u.tr() * meas(G), u * uex * meas(G) ); + solver.compute(A.matrix()); + solVector_ini = solVector = solver.solve(A.rhs()); + + gsMatrix<> fluxData(2,fluxControlPoints.transpose().cols()), tmp2; + // Write initial data + if (participant.requiresWritingCheckpoint()) + { + for (index_t k=0; k!=fluxControlPoints.rows(); k++) + { + gsWarn<<"Write the flux here!!!\n"; + tmp2 = ev.eval( - jac(u_sol) * nv(G).normalized(),fluxControlPoints.transpose().col(k)); + gsInfo << "Got here\n"; + fluxData(0,k) = tmp2.at(0); + } + gsDebugVar(fluxData); + // gsDebugVar(result); + participant.writeData(FluxControlPointMesh, FluxControlPointData, fluxControlPointsIDs, fluxData); + + } + participant.readData(GeometryControlPointMesh,TemperatureData,geometryControlPointIDs,tempData); + tempMesh.patch(0).coefs()=tempData.row(0).transpose(); + gsDebugVar(tempMesh.patch(0).coefs()); + g_C = &tempMesh.patch(0); //Update the boundary condition for the east coupled boundary + A.initSystem(); + + // Assemble stiffness matrix $K = \int_{\Omega} \kappa \nabla \varphi_i \cdot \nabla \varphi_j , d\Omega$ + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * uex * meas(G) ); + gsSparseMatrix<> K = A.matrix(); + // Discretize time, the system becomes: $(M + \Delta t K) u^{n+1} = M u^n + \Delta t f$ + + gsVector<> F = dt*A.rhs() + M*solVector; + gsVector<> F0 = F; + gsVector<> F_checkpoint = F; + gsMatrix<> solVector_checkpoint = solVector; + + gsParaviewCollection collection("solution_dirichlet"); + collection.options().setSwitch("plotElements",true); + gsParaviewCollection exact_collection("exact_solution_dirichlet"); + + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + real_t t_checkpoint = 0; + if (plot) + { + std::string fileName = "solution_dirichlet" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_dirichlet" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + // fileName = "exact_solution_dirichlet" + util::to_string(timestep); + // ev.writeParaview( uex , G, fileName); + // for (size_t p=0; p!=patches.nPatches(); p++) + // { + // fileName = "exact_solution_dirichlet" + util::to_string(timestep) + std::to_string(p); + // exact_collection.addTimestep(fileName,time,".vts"); + // } + + // ev.writeParaview( u_sol , G, "initial_solution_dirichlet"); + + } + while (participant.isCouplingOngoing()) + { + + u_ex = gsFunctionExpr<>("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); + u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * ff * meas(G) ); + K = A.matrix(); + F = dt*A.rhs() + M*solVector; + // save checkpoint + if (participant.requiresWritingCheckpoint()) + { + gsDebugVar("Write checkpoint"); + F_checkpoint = F0; + t_checkpoint = time; + timestep_checkpoint = timestep; + solVector_checkpoint = solVector; + } + + // potentially adjust non-matching timestep sizes + dt = std::min(dt,precice_dt); + + // solve gismo timestep + gsInfo << "Solving timestep " << timestep*dt << "..."; + solVector = solver.compute(M + dt*K).solve(F); + gsInfo<<"Finished\n"; + // write heat fluxes to interface + // gsMatrix<> result(2,fluxControlPoints.transpose().cols()), tmp; + + for (index_t k=0; k!=fluxControlPoints.transpose().cols(); k++) + { + gsWarn<<"Write the flux here!!!\n"; + //Calculate the heat flux $q=-k\nabla u$, here we use jac(u_sol) to calculate $\nabla u$ + // Use normal vector to calculate the flux in the normal direction + tmp2 = ev.eval(-jac(u_sol) * nv(G).normalized(),fluxControlPoints.transpose().col(k)); + gsInfo << "Got here\n"; + fluxData(0,k) = tmp2.at(0); + } + participant.writeData(FluxControlPointMesh, FluxControlPointData, fluxControlPointsIDs, fluxData); + participant.readData(GeometryControlPointMesh,TemperatureData,geometryControlPointIDs,tempData); + gsDebugVar(tempData); + tempMesh.patch(0).coefs()=tempData.row(0).transpose(); + gsDebugVar(tempData.row(0).transpose()); + gsDebugVar(tempMesh.patch(0).coefs()); + g_C = &tempMesh.patch(0); + // g_C = &tempMesh.patch(0); //Update the boundary condition for the east coupled boundary + + // do the coupling + precice_dt = participant.advance(dt); + + // advance variables + time += dt; + timestep += 1; + F0 = F; + + if (participant.requiresReadingCheckpoint()) + { + gsDebugVar("Read checkpoint"); + F0 = F_checkpoint; + time = t_checkpoint; + timestep = timestep_checkpoint; + solVector = solVector_checkpoint; + } + else + { + if (timestep % plotmod==0 && plot) + { + std::string fileName = "solution_dirichlet" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_dirichlet" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + + } + } + } + + if (plot) + { + collection.save(); + exact_collection.save(); + } + return EXIT_SUCCESS; +} diff --git a/examples/partitioned-heat-conduction-IGA-IGA/neumann/partitioned-heat-conduction-IGA-Neumann.cpp b/examples/partitioned-heat-conduction-IGA-IGA/neumann/partitioned-heat-conduction-IGA-Neumann.cpp new file mode 100644 index 0000000..ebd7ded --- /dev/null +++ b/examples/partitioned-heat-conduction-IGA-IGA/neumann/partitioned-heat-conduction-IGA-Neumann.cpp @@ -0,0 +1,385 @@ +/// This is the fluid-structure interaction benchmark FSI2 from this paper: +/// "Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow" +/// Stefan Turek and Jaroslav Hron, , 2006. +/// Author(s): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) +/// +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +// #include + +using namespace gismo; + + +int main(int argc, char* argv[]) +{ + gsInfo << "Testing the two-way fluid-structure interaction solver in 2D.\n"; + + //=====================================// + // Input // + //=====================================// + + // beam parameters + bool plot = true; + index_t plotmod = 1; + index_t numUniRef = 3; + index_t numElevate = 0; + short_t side = 0; + real_t alpha = 3; + real_t beta = 1.2; + real_t time = 0; + real_t k_temp = 1; + + + std::string precice_config("../precice_config.xml"); + + // minimalistic user interface for terminal + gsCmdLine cmd("Testing the two-way fluid-structure interaction solver in 2D."); + cmd.addInt("r","refine","Number of uniform refinement applications",numUniRef); + // cmd.addReal("t","time","Time span, sec",timeSpan); + // cmd.addReal("s","step","Time step",timeStep); + // cmd.addInt("p","points","Number of points to plot to Paraview",numPlotPoints); + // cmd.addInt("v","verbosity","Amount of info printed to the prompt: 0 - none, 1 - crucial, 2 - all",verbosity); + cmd.addString( "c", "config", "PreCICE config file", precice_config ); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + GISMO_ASSERT(gsFileManager::fileExists(precice_config),"No precice config file has been defined"); + + //=============================================// + // Scanning geometry and creating bases // + //=============================================// + + // scanning geometry + gsMultiPatch<> patches; + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(1.0,0.0,2.0,1.0)); + + // creating bases + for (index_t i = 0; i < numUniRef; ++i) + patches.uniformRefine(); + gsMultiBasis<> basisTemperature(patches); + + // Set external heat-flux to zero + gsConstantFunction<> f(beta-2-2*alpha,2); + gsFunctionExpr<> u_ex("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); + + + //=============================================// + // Define preCICE setup // + //=============================================// + + std::string participantName = "Neumann"; + gsPreCICE participant(participantName, precice_config); + + // Mesh provided by solid + std::string GeometryKnotMesh = "Geometry-Knot-Mesh"; + std::string GeometryControlPointMesh = "Geometry-Control-Point-Mesh"; + + // Mesh provided by Dirichlet + std::string FluxKnotMesh = "Flux-Knot-Mesh"; + std::string FluxControlPointMesh = "Flux-Control-Point-Mesh"; + + // Data provided by Neumann + std::string TemperatureData = "Temperature-Data"; + + // Data provided by Dirichlet + std::string FluxControlPointData = "Flux-Control-Point-Data"; + + // Setup bounding box onto the force mesh + gsMatrix<> bbox(2,2); + bbox << -1e300, 1e300, // X dimension limits + -1e300, 1e300; // Y dimension limits + + bbox.transposeInPlace(); + bbox.resize(1,bbox.rows()*bbox.cols()); + participant.setMeshAccessRegion(FluxControlPointMesh, bbox); + + // Setup geometry control point mesh + gsVector geometryControlPointIDs; + gsMatrix<> geometryControlPoints = patches.patch(0).coefs().transpose(); + + std::vector couplingInterface(1); + couplingInterface[0] = patchSide(0,boundary::south); + // std::vector::uPtr> boundaries(couplingInterface.size()); + + gsMultiBasis<> temperatureBases; + std::vector::uPtr> couplingBoundaries(couplingInterface.size()); + couplingBoundaries[0] = patches.patch(0).boundary(couplingInterface[0].side()); + + gsMatrix<> temperatureControlPoints = couplingBoundaries[0]->coefs(); + participant.addMesh(GeometryControlPointMesh, temperatureControlPoints.transpose(), geometryControlPointIDs); + + // gsMatrix<> boundaryPoints(temperatureControlPoints.cols(),temperatureControlPoints.rows()); + // boundaryPoints.row(0) = temperatureControlPoints.col(0).transpose(); + + // gsDebugVar(boundaryPoints); + gsMatrix<> boundaryPoints = temperatureControlPoints; + boundaryPoints.col(0).array() -= 1; + gsDebugVar(boundaryPoints); + + gsVector<> temperatureKnotMatrix = knotsToVector(couplingBoundaries[0]->basis()); + gsDebugVar(temperatureKnotMatrix); + participant.addMesh(GeometryKnotMesh, temperatureKnotMatrix.transpose()); + + + + // Setup the geometry knot mesh + + //participant.initialize(); + + // Initialize preCICE (send solid mesh to API) + real_t precice_dt = participant.initialize(); + + + //Get the force mesh from direct-access="true"direct-access="true"the API + gsVector fluxKnotIDs; + gsMatrix<> fluxKnots; + participant.getMeshVertexIDsAndCoordinates(FluxKnotMesh,fluxKnotIDs,fluxKnots); + + gsDebugVar(fluxKnots); + + // Gives a full tensor product basis + + gsBasis<> * basis = knotMatrixToBasis(fluxKnots.row(0)).get(); + + gsDebugVar(basis->size()); + + + gsMatrix<> fluxControlPoints(basis->size(),1); + fluxControlPoints.setZero(); + + + // Gives the coefficients of the control points ONLY ON THE BOUNDARIES + + // IMPORTANT: THE control points now are all control points!!! + gsVector fluxControlPointIDs; + participant.getMeshVertexIDsAndCoordinates(FluxControlPointMesh, fluxControlPointIDs,fluxControlPoints); + gsDebugVar(fluxControlPoints); + + + + + // // Step 2: Regenerate the geometry + gsMultiPatch<> fluxMesh; //Geometry object belongs to gsFunctionSet + fluxMesh.addPatch(give(basis->makeGeometry(fluxControlPoints.transpose()))); + + gsDebugVar(fluxMesh.patch(0).coefs()); + + + // NOTE: forceMesh should have domainDim 2!! + //=============================================// + // Setting loads and boundary conditions // + //=============================================// + // Setup the boundary condition for the Neumann side + gsBoundaryConditions<> bcInfo; + gsConstantFunction<> g_D(0,2); // Dirichlet + bcInfo.addCondition(0, boundary::south, condition_type::neumann , &fluxMesh.patch(0)); + bcInfo.addCondition(0, boundary::east, condition_type::dirichlet , &u_ex); + bcInfo.addCondition(0, boundary::west, condition_type::dirichlet , &u_ex); + bcInfo.addCondition(0, boundary::north, condition_type::dirichlet , &u_ex); + + bcInfo.setGeoMap(patches); + gsDebugVar(bcInfo); + + // ---------------------------------------------------------------------------------------------- + + // Generate system matrix and load vector + // gsInfo << "Assembling mass and stiffness...\n"; + // + // ---------------------------------------------------------------------------------------------- + typedef gsExprAssembler<>::geometryMap geometryMap; + typedef gsExprAssembler<>::space space; + typedef gsExprAssembler<>::solution solution; + + gsExprAssembler<> A(1,1); + + gsInfo<<"Active options:\n"<< A.options() <<"\n"; + + A.setIntegrationElements(basisTemperature); + + gsExprEvaluator<> ev(A); + + + // Set the geometry map + geometryMap G = A.getMap(patches); + + // Set the discretization space + space u = A.getSpace(basisTemperature); + + // Set the source term + auto ff = A.getCoeff(f, G); + + // Set the solution + gsMatrix<> solVector, solVector_ini; + solution u_sol = A.getSolution(u, solVector); + solution u_ini = A.getSolution(u, solVector); + + u.setup(bcInfo, dirichlet::homogeneous, 0); // Add boundary condition + A.initSystem(); + gsDebugVar(A.numDofs()); + A.assemble( u * u.tr() * meas(G)); + gsSparseMatrix<> M = A.matrix(); + + // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner + gsSparseSolver<>::CGDiagonal solver; + + real_t dt = 0.1; + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + auto uex = A.getCoeff(u_ex, G); + + // RHS of the projection + u.setup(bcInfo, dirichlet::l2Projection, 0); + A.initSystem(); + A.assemble( u * u.tr() * meas(G), u * uex * meas(G) ); + solver.compute(A.matrix()); + solVector_ini = solVector = solver.solve(A.rhs()); + + gsMatrix<> results(2,boundaryPoints.transpose().cols()), tmp2; + + + // Write initial data + if (participant.requiresWritingCheckpoint()) + { + + for (index_t k=0; k!=boundaryPoints.transpose().cols(); k++) + { + gsInfo << "Write Temperature Here! \n"; + tmp2 = ev.eval(u_sol,boundaryPoints.transpose().col(k)); + gsDebugVar(tmp2); + results(0,k) = tmp2.at(0); + } + gsDebugVar(results); + participant.writeData(GeometryControlPointMesh, TemperatureData, geometryControlPointIDs, results); + } + participant.readData(FluxControlPointMesh,FluxControlPointData,fluxControlPointIDs,fluxControlPoints); + fluxMesh.patch(0).coefs() = fluxControlPoints; + gsDebugVar(fluxControlPoints); + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * uex * meas(G) ); + gsSparseMatrix<> K = A.matrix(); + + gsVector<> F = dt*A.rhs() + M*solVector; + gsVector<> F0 = F; + gsVector<> F_checkpoint = F; + gsMatrix<> solVector_checkpoint = solVector; + + gsParaviewCollection collection("solution_neumann"); + collection.options().setSwitch("plotElements",true); + gsParaviewCollection exact_collection("exact_solution_neumann"); + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + real_t t_checkpoint = 0; + if (plot) + { + std::string fileName = "solution_neumann" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_neumann" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + ev.writeParaview( u_sol , G, "initial_solution_neumann"); + + } + + while (participant.isCouplingOngoing()) + { + + u_ex = gsFunctionExpr<>("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); + u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: + + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * ff * meas(G) ); + K = A.matrix(); + // auto g_Neumann = A.getBdrFunction(G); + // A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); + F = dt*A.rhs() + M*solVector; + // save checkpoint + if (participant.requiresWritingCheckpoint()) + { + gsDebugVar("Write checkpoint"); + F_checkpoint = F0; + t_checkpoint = time; + timestep_checkpoint = timestep; + solVector_checkpoint = solVector; + } + + // potentially adjust non-matching timestep sizes + dt = std::min(dt,precice_dt); + + // solve gismo timestep + gsInfo << "Solving timestep " << timestep*dt << "..."; + solVector = solver.compute(M + dt*K).solve(F); + gsInfo<<"Finished\n"; + // write heat fluxes to interface + // gsMatrix<> result(2,fluxControlPoints.transpose().cols()), tmp; + // + // + gsDebugVar(fluxControlPoints); + for (index_t k=0; k!=boundaryPoints.transpose().cols(); k++) + { + gsInfo << "Write Temperature Here! \n"; + tmp2 = ev.eval(u_sol,boundaryPoints.transpose().col(k).transpose()); + results(0,k) = tmp2.at(0); + } + gsDebugVar(results); + participant.writeData(GeometryControlPointMesh, TemperatureData, geometryControlPointIDs, results); + participant.readData(FluxControlPointMesh,FluxControlPointData,fluxControlPointIDs,fluxControlPoints); + fluxMesh.patch(0).coefs() = fluxControlPoints; + gsDebugVar(fluxControlPoints); + // do the coupling + precice_dt = participant.advance(dt); + + // advance variables + time += dt; + timestep += 1; + F0 = F; + + if (participant.requiresReadingCheckpoint()) + { + gsDebugVar("Read checkpoint"); + F0 = F_checkpoint; + time = t_checkpoint; + timestep = timestep_checkpoint; + solVector = solVector_checkpoint; + } + else + { + if (timestep % plotmod==0 && plot) + { + std::string fileName = "solution_neumann" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_neumann" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + + } + } + } + + if (plot) + { + collection.save(); + exact_collection.save(); + } + + return EXIT_SUCCESS; +} diff --git a/examples/partitioned-heat-conduction-IGA-IGA/precice-config.xml b/examples/partitioned-heat-conduction-IGA-IGA/precice-config.xml new file mode 100644 index 0000000..899d1b4 --- /dev/null +++ b/examples/partitioned-heat-conduction-IGA-IGA/precice-config.xml @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/partitioned-heat-conduction-direct/.DS_Store b/examples/partitioned-heat-conduction-direct/.DS_Store new file mode 100644 index 0000000..03d268d Binary files /dev/null and b/examples/partitioned-heat-conduction-direct/.DS_Store differ diff --git a/examples/partitioned-heat-conduction/gismo-left/clean.sh b/examples/partitioned-heat-conduction-direct/gismo-left/clean.sh similarity index 100% rename from examples/partitioned-heat-conduction/gismo-left/clean.sh rename to examples/partitioned-heat-conduction-direct/gismo-left/clean.sh diff --git a/examples/partitioned-heat-conduction/gismo-left/run.sh b/examples/partitioned-heat-conduction-direct/gismo-left/run.sh similarity index 100% rename from examples/partitioned-heat-conduction/gismo-left/run.sh rename to examples/partitioned-heat-conduction-direct/gismo-left/run.sh diff --git a/examples/partitioned-heat-conduction/gismo-right/clean.sh b/examples/partitioned-heat-conduction-direct/gismo-right/clean.sh similarity index 100% rename from examples/partitioned-heat-conduction/gismo-right/clean.sh rename to examples/partitioned-heat-conduction-direct/gismo-right/clean.sh diff --git a/examples/partitioned-heat-conduction/gismo-right/run.sh b/examples/partitioned-heat-conduction-direct/gismo-right/run.sh similarity index 100% rename from examples/partitioned-heat-conduction/gismo-right/run.sh rename to examples/partitioned-heat-conduction-direct/gismo-right/run.sh diff --git a/examples/partitioned-heat-conduction-direct/partitioned-heat-conduction-direct.cpp b/examples/partitioned-heat-conduction-direct/partitioned-heat-conduction-direct.cpp new file mode 100644 index 0000000..776ae84 --- /dev/null +++ b/examples/partitioned-heat-conduction-direct/partitioned-heat-conduction-direct.cpp @@ -0,0 +1,407 @@ +/** @file heat-equation-coupling.cpp + + @brief Heat equation participant for a double coupled heat equation + + 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): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) +*/ + +#include +#include +#include + +using namespace gismo; + +int main(int argc, char *argv[]) +{ + //! [Parse command line] + bool plot = false; + index_t plotmod = 1; + index_t numRefine = 5; + index_t numElevate = 0; + short_t side = 0; + std::string precice_config("../precice_config.xml"); + + gsCmdLine cmd("Coupled heat equation using 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("s","side", "Patchside of interface", side); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + //! [Read input file] + + gsMultiPatch<> patches; + if (side==0) //left + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,0.0,1.0,1.0)); + else if (side==1) //right + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(1.0,0.0,2.0,1.0)); + else + GISMO_ERROR("Side unknown"); + + real_t alpha = 3; + real_t beta = 1.3; + real_t time = 0; + real_t k_temp = 1; + + // Set external heat-flux to zero + gsConstantFunction<> f(beta-2-2*alpha,2); + gsFunctionExpr<> u_ex("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); + + gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) + + 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"; + +// ---------------------------------------------------------------------------------------------- + typedef gsExprAssembler<>::geometryMap geometryMap; + typedef gsExprAssembler<>::space space; + typedef gsExprAssembler<>::solution solution; + + gsExprAssembler<> A(1,1); + + gsInfo<<"Active options:\n"<< A.options() <<"\n"; + + A.setIntegrationElements(bases); + + gsExprEvaluator<> ev(A); + +// ---------------------------------------------------------------------------------------------- + boxSide couplingSide; + if (side==0) //left + couplingSide = boxSide(2); + else if (side==1) //right + couplingSide = boxSide(1); + else + GISMO_ERROR("Side unknown"); + + patchSide couplingInterface(0,couplingSide); + typename gsBasis::domainIter domIt = bases.basis(couplingInterface.patch).makeDomainIterator(couplingInterface.side()); + index_t rows = patches.targetDim(); + gsMatrix<> nodes; + // Start iteration over elements + gsVector<> tmp; + index_t k=0; + + gsOptionList quadOptions = A.options(); + // NEED A DIFFERENT RULE FOR dirichlet::interpolation --> see: gsGaussRule bdQuRule(basis, 1.0, 1, iter->side().direction()); + /* + quadOptions.addInt("quRule","Quadrature rule [1:GaussLegendre,2:GaussLobatto]",1); + quadOptions.addReal("quA","Number of quadrature points: quA*deg + quB",1.0); + quadOptions.addInt("quB","Number of quadrature points: quA*deg + quB",1); + */ + + index_t quadSize = 0; + typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT + for (; domIt->good(); domIt->next(), k++ ) + { + QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); + quadSize+=QuRule->numNodes(); + } + gsMatrix<> uv(rows,quadSize); + gsMatrix<> xy(rows,quadSize); + + index_t offset = 0; + + for (domIt->reset(); domIt->good(); domIt->next(), k++ ) + { + QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); + // Map the Quadrature rule to the element + QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(), + nodes, tmp); + uv.block(0,offset,rows,QuRule->numNodes()) = nodes; + + gsMatrix<> tmp2; + patches.patch(couplingInterface.patch).eval_into(nodes,tmp2); + xy.block(0,offset,rows,QuRule->numNodes()) = patches.patch(couplingInterface.patch).eval(nodes); + offset += QuRule->numNodes(); + } + + GISMO_ASSERT(side==0 || side==1,"Side must be west or east"); + std::string readName = (side == 0) ? "Dirichlet" : "Neumann"; + std::string writeName = (side == 0) ? "Neumann" : "Dirichlet"; + + gsPreCICE interface(membername, precice_config); + std::string readMeshName = readName + "-Mesh"; + std::string writeMeshName = writeName + "-Mesh"; + interface.addMesh(readMeshName,xy); + + interface.requiresInitialData(); + + real_t precice_dt = interface.initialize(); + + std::string tempName = "Temperature"; + std::string fluxName = "Heat-Flux"; + + gsMatrix<> bbox(2,2); + bbox<<0.9,1.1,-0.1,1.1; + interface.setMeshAccessRegion(writeMeshName,bbox); + + std::vector writeIDs; + gsMatrix<> writePoints; + interface.getMeshVertexIDsAndCoordinates(writeMeshName,writeIDs,writePoints); + +/* + TO DO: + * Write the evaluation on the other coordinates on the other IDs +*/ +// ---------------------------------------------------------------------------------------------- + + gsBoundaryConditions<> bcInfo; + gsPreCICEFunction g_CD(&interface,readMeshName,(side==0 ? tempName : fluxName),patches,1); + gsPreCICEFunction g_CN(&interface,readMeshName,(side==0 ? tempName : fluxName),patches,1); + gsFunction<> * g_C = &u_ex; + + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &u_ex, 0, false, 0); + bcInfo.addCondition(0, boundary::north, condition_type::dirichlet, &u_ex, 0, false, 0); + if (side==0) + { + bcInfo.addCondition(0, couplingSide,condition_type::dirichlet , g_C, 0, false, 0); + bcInfo.addCondition(0, couplingSide.opposite(),condition_type::dirichlet, &u_ex, 0, false, 0); + } + else + { + bcInfo.addCondition(0, couplingSide,condition_type::neumann , &g_CN); + bcInfo.addCondition(0, couplingSide.opposite(),condition_type::dirichlet, &u_ex, 0, false, 0); + } + + bcInfo.setGeoMap(patches); + // gsDebugVar(bcInfo); + +// ---------------------------------------------------------------------------------------------- + + // Generate system matrix and load vector + // gsInfo << "Assembling mass and stiffness...\n"; + + // Set the geometry map + geometryMap G = A.getMap(patches); + + // Set the discretization space + space u = A.getSpace(bases); + + // Set the source term + auto ff = A.getCoeff(f, G); + + // Set the solution + gsMatrix<> solVector, solVector_ini; + solution u_sol = A.getSolution(u, solVector); + solution u_ini = A.getSolution(u, solVector); + + u.setup(bcInfo, dirichlet::homogeneous, 0); + A.initSystem(); + gsDebugVar(A.numDofs()); + A.assemble( u * u.tr() * meas(G)); + gsSparseMatrix<> M = A.matrix(); + + + // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner + gsSparseSolver<>::CGDiagonal solver; + + real_t dt = 0.01; + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + auto uex = A.getCoeff(u_ex, G); + // RHS of the projection + u.setup(bcInfo, dirichlet::l2Projection, 0); + A.initSystem(); + A.assemble( u * u.tr() * meas(G), u * uex * meas(G) ); + solver.compute(A.matrix()); + solVector_ini = solVector = solver.solve(A.rhs()); + + gsMatrix<> result(1,uv.cols()), tmp2; + // Write initial data + if (interface.requiresWritingCheckpoint()) + { + for (index_t k=0; k!=uv.cols(); k++) + { + if (side==0) + { + gsWarn<<"Write the flux here!!!\n"; + tmp2 = ev.eval( - jac(u_sol) * nv(G).normalized(),writePoints.col(k)); + result(0,k) = tmp2.at(0); + } + else + { + tmp2 = ev.eval(u_sol,writePoints.col(k)); + result(0,k) = tmp2.at(0); + } + } + // gsDebugVar(result); + interface.writeData(writeMeshName,writeName=="Dirichlet" ? fluxName : tempName,writeIDs,result); + } + // interface.initialize_data(); + + // interface.readBlockScalarData(readMeshName,side==0 ? tempName : fluxName,xy,result); + // gsDebugVar(result); + + // Initialize the RHS for assembly + if (side==0) + g_C = &g_CD; + gsDebugVar(bcInfo); + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * uex * meas(G) ); + gsSparseMatrix<> K = A.matrix(); + + // Assemble the RHS + gsVector<> F = dt*A.rhs() + M*solVector; + gsVector<> F0 = F; + gsVector<> F_checkpoint = F; + gsMatrix<> solVector_checkpoint = solVector; + + gsParaviewCollection collection("solution"); + gsParaviewCollection exact_collection("exact_solution"); + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + real_t t_checkpoint = 0; + if (plot) + { + std::string fileName = "solution_" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + + // fileName = "exact_solution_" + util::to_string(timestep); + // ev.writeParaview( uex , G, fileName); + // for (size_t p=0; p!=patches.nPatches(); p++) + // { + // fileName = "exact_solution_" + util::to_string(timestep) + std::to_string(p); + // exact_collection.addTimestep(fileName,time,".vts"); + // } + + ev.writeParaview( u_sol , G, "initial_solution"); + + } + + // time += precice_dt; + while (interface.isCouplingOngoing()) + { + // read temperature from interface + // if (interface.isReadDataAvailable()) + // { + u_ex = gsFunctionExpr<>("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); + u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: + + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * ff * meas(G) ); + K = A.matrix(); + auto g_Neumann = A.getBdrFunction(G); + A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); + F = dt*A.rhs() + M*solVector; + + // gsMatrix<> result; + // interface.readBlockScalarData(readMeshName,side==0 ? tempName : fluxName,xy,result); + // gsDebugVar(result); + // } + + // save checkpoint + if (interface.requiresWritingCheckpoint()) + { + gsDebugVar("Write checkpoint"); + F_checkpoint = F0; + t_checkpoint = time; + timestep_checkpoint = timestep; + solVector_checkpoint = solVector; + } + + // potentially adjust non-matching timestep sizes + dt = std::min(dt,precice_dt); + + // solve gismo timestep + gsInfo << "Solving timestep " << timestep*dt << "..."; + solVector = solver.compute(M + dt*K).solve(F); + gsInfo<<"Finished\n"; + // write heat fluxes to interface + gsMatrix<> result(1,uv.cols()), tmp; + for (index_t k=0; k!=uv.cols(); k++) + { + // gsDebugVar(ev.eval(nv(G),uv.col(k))); + // tmp = ev.eval(k_temp * igrad(u_sol,G),uv.col(k)); + // Only exchange y component + // result(0,k) = -tmp.at(1); + // + // + if (side==0) + { + gsWarn<<"Write the flux here!!!\n"; + tmp = ev.eval(jac(u_sol) * nv(G).normalized(),uv.col(k)); + result(0,k) = tmp.at(0); + } + else + { + tmp = ev.eval(u_sol,uv.col(k)); + result(0,k) = tmp.at(0); + } + } + // TODO + interface.writeData(readMeshName,side==0 ? fluxName : tempName,xy,result); + + // do the coupling + precice_dt = interface.advance(dt); + + // advance variables + time += dt; + timestep += 1; + F0 = F; + + if (interface.requiresReadingCheckpoint()) + { + gsDebugVar("Read checkpoint"); + F0 = F_checkpoint; + time = t_checkpoint; + timestep = timestep_checkpoint; + solVector = solVector_checkpoint; + } + else + { + if (timestep % plotmod==0 && plot) + { + std::string fileName = "solution_" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + + // fileName = "exact_solution_" + util::to_string(timestep); + // ev.writeParaview( uex , G, fileName); + // for (size_t p=0; p!=patches.nPatches(); p++) + // { + // fileName = "exact_solution_" + util::to_string(timestep) + std::to_string(p); + // exact_collection.addTimestep(fileName,time,".vts"); + // } + } + } + } + + if (plot) + { + collection.save(); + exact_collection.save(); + } + + + return EXIT_SUCCESS; +} diff --git a/examples/partitioned-heat-conduction-direct/precice-config.xml b/examples/partitioned-heat-conduction-direct/precice-config.xml new file mode 100644 index 0000000..46d6bdf --- /dev/null +++ b/examples/partitioned-heat-conduction-direct/precice-config.xml @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/partitioned-heat-conduction-vertex-vertex/.DS_Store b/examples/partitioned-heat-conduction-vertex-vertex/.DS_Store new file mode 100644 index 0000000..d848885 Binary files /dev/null and b/examples/partitioned-heat-conduction-vertex-vertex/.DS_Store differ diff --git a/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/clean.sh b/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/clean.sh new file mode 100755 index 0000000..708528e --- /dev/null +++ b/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_gismo . diff --git a/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/run.sh b/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/run.sh new file mode 100755 index 0000000..7630e55 --- /dev/null +++ b/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/run.sh @@ -0,0 +1,4 @@ +#!/bin/bash +set -e -u + +../../../../../build/bin/partitioned-heat-conduction -c ../precice-config.xml --plot -s 0 diff --git a/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/clean.sh b/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/clean.sh new file mode 100755 index 0000000..708528e --- /dev/null +++ b/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_gismo . diff --git a/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/run.sh b/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/run.sh new file mode 100755 index 0000000..026edd8 --- /dev/null +++ b/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/run.sh @@ -0,0 +1,4 @@ +#!/bin/bash +set -e -u + +../../../../../build/bin/partitioned-heat-conduction -c ../precice-config.xml --plot -s 1 diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction.cpp b/examples/partitioned-heat-conduction-vertex-vertex/partitioned-heat-conduction.cpp similarity index 100% rename from examples/partitioned-heat-conduction/partitioned-heat-conduction.cpp rename to examples/partitioned-heat-conduction-vertex-vertex/partitioned-heat-conduction.cpp diff --git a/examples/partitioned-heat-conduction/precice-config.xml b/examples/partitioned-heat-conduction-vertex-vertex/precice-config.xml similarity index 100% rename from examples/partitioned-heat-conduction/precice-config.xml rename to examples/partitioned-heat-conduction-vertex-vertex/precice-config.xml diff --git a/examples/perpendicular-flap-IGA-IGA/perpendicular-flap-IGA-IGA-fluid.cpp b/examples/perpendicular-flap-IGA-IGA/perpendicular-flap-IGA-IGA-fluid.cpp new file mode 100644 index 0000000..67b2ab2 --- /dev/null +++ b/examples/perpendicular-flap-IGA-IGA/perpendicular-flap-IGA-IGA-fluid.cpp @@ -0,0 +1,282 @@ +/** @file test-communication-IGA-IGA-fluid.cpp + * + * @brief test communicate spline through preCICE. To see if the spline force function can be reconstructed + * + * This file is a part of 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 (TU Delft, 2023-...) + * + **/ + +#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; + bool write = false; + bool get_readTime = false; + bool get_writeTime = false; + index_t plotmod = 1; + index_t loadCase = 0; + std::string precice_config("../precice_config.xml"); + + gsCmdLine cmd("Coupled heat equation using PreCICE."); + cmd.addString( "c", "config", "PreCICE config file", precice_config ); + cmd.addInt("l","loadCase", "Load case: 0=constant load, 1='spring' load", loadCase); + 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); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + GISMO_ASSERT(gsFileManager::fileExists(precice_config),"No precice config file has been defined"); + //! [Read input file] + std::string participantName = "Fluid"; + gsPreCICE participant(participantName, precice_config); + + + /* + * Data initialization + * + * This participant receives mesh information (knot vector, control points) from the solid, + * and it creates its own mesh based on the information. + * And writes pressure, reads displacements from the solid participant. + * The follow meshes and data are made available: + * + * - Meshes: + * + GeometryControlPointMesh: This mesh contains the control points as mesh vertices. + * + GeometryKnotMesh: This mesh contains the knots as mesh vertices. + * + ForceKnotMesh: This mesh contains the knots as mesh vertices for force mesh. + * + ForceControlPointMesh: This mesh contains the control points as mesh vertices for force mesh. + * + * + * - Data: + * + GeometryControlPointData: This data is defined on the ControlPointMesh and stores the displacement of the control points + * + ForceControlPointData: This data is defined on the ForceControlPointMesh and stores the change of pressure/stress + */ + + + // Meshes + std::string GeometryControlPointMesh = "GeometryControlPointMesh"; + std::string GeometryKnotMesh = "GeometryKnotMesh"; + + + std::string ForceKnotMesh = "ForceKnotMesh"; + std::string ForceControlPointMesh = "ForceControlPointMesh"; + + + // Data + std::string GeometryControlPointData = "GeometryControlPointData"; + std::string ForceControlPointData = "ForceControlPointData"; + // std::string GeometryKnotData = "GeometryKnotData"; + // std::string ForceKnotData = "ForceKnotMeshData"; + + gsMatrix<> bbox(3,2); + bbox << -1e300, 1e300, // X dimension limits + -1e300, 1e300, // Y dimension limits + -1e300, 1e300; // Z dimension limits + bbox.transposeInPlace(); + bbox.resize(1,bbox.rows()*bbox.cols()); + + // Step 2: Add force mesh + + // Step 2.1: Define force knot mesh + gsVector forceKnotIDs; + gsKnotVector<> kv(0,1,10,3,1); //begin, end, # interiors, mult, by default =1 + gsTensorBSplineBasis<2, real_t> forceBasis(kv,kv); + gsMatrix<> forceControlPoints(forceBasis.size(), 3); + forceControlPoints.setZero(); + forceControlPoints.col(2).setConstant(-5e3); + forceControlPoints.transposeInPlace(); + + gsMatrix<> forceKnots = knotsToMatrix(forceBasis); + + + + + // gsMultiBasis<> bases(forceMesh); + // gsMatrix<> forceKnots = knotsToMatrix(bases.basis(0)); + // gsDebugVar(forceKnots.dim()); + // + + + + // Regenerate the geometry in another solver + // forceMesh.addPatch(give(forceBasis.makeGeometry(forceControlPoints))); + + // Step 2.1: Define force control point mesh + gsVector forceControlPointIDs; + participant.addMesh(ForceKnotMesh,forceKnots,forceKnotIDs); + participant.addMesh(ForceControlPointMesh,forceControlPoints,forceControlPointIDs); + participant.setMeshAccessRegion(GeometryControlPointMesh, bbox); + + real_t precice_dt = participant.initialize(); + + // Step 1: Collect geometry from PreCICE + gsVector geometryKnotIDs; + gsMatrix<> geometryKnots; + participant.getMeshVertexIDsAndCoordinates(GeometryKnotMesh, geometryKnotIDs, geometryKnots); + gsDebugVar(geometryKnots); + + + gsBasis<> * geometryKnotBasis = knotMatrixToBasis(geometryKnots).get(); + + gsVector geometryControlPointIDs; + gsMatrix<> geometryControlPoints; + participant.getMeshVertexIDsAndCoordinates(GeometryControlPointMesh, geometryControlPointIDs, geometryControlPoints); + gsDebugVar(geometryControlPoints.dim()); + + gsMultiPatch<> mp, deformation; + gsDebugVar(*geometryKnotBasis); + gsDebugVar(geometryKnotBasis->size()); + gsDebugVar(geometryControlPoints.rows()); + mp.addPatch(give(geometryKnotBasis->makeGeometry(geometryControlPoints.transpose()))); + + deformation = mp; + + + real_t t=0; + real_t dt = precice_dt; + real_t t_read = 0; + real_t t_write = 0; + + index_t timestep = 0; + + // Define the solution collection for Paraview + gsParaviewCollection collection("./output/solution"); + + gsMatrix<> points(2,1); + points.col(0)<<0.5,1; + + gsStructuralAnalysisOutput writer("./output/pointData.csv",points); + writer.init({"x","y","z"},{"time"}); + + // Time integration loop + while(participant.isCouplingOngoing()) + { + + // Read control point displacements + participant.readData(GeometryControlPointMesh, GeometryControlPointData, geometryControlPointIDs, geometryControlPoints); + deformation.patch(0).coefs() = geometryControlPoints.transpose(); + if (get_readTime) + t_read += participant.readTime(); + + /* + * Two projection approaches: + * - gsQuasiInterpolate::localIntpl(forceBasis, deformation.patch(0), coefs); + * - gsL2Projection::projectFunction(forceBasis, deformation, mp, coefs) + * + * Check if forceBasis + coefs gives the same as deformation.patch(0) + */ + if (loadCase == 0) + { + gsInfo << "Load case 0: Constant load\n"; + } + else if (loadCase == 1) + { + gsInfo << "Load case 1: Spring load\n"; + if(timestep > 0) + { + gsMatrix<> coefs; + gsL2Projection::projectFunction(forceBasis, deformation, mp, coefs); + coefs.resize(forceControlPoints.rows(),forceControlPoints.cols()); + forceControlPoints.row(2) = -coefs.row(2); + } + } + else + { + GISMO_ERROR("Load case "< pointDataMatrix; + gsMatrix<> otherDataMatrix(1,1); + if (participant.requiresWritingCheckpoint()) + { + gsInfo<<"Reading Checkpoint\n"; + } + + // do the coupling + precice_dt = participant.advance(dt); + + dt = std::min(precice_dt, dt); + + if (participant.requiresReadingCheckpoint()) + { + gsInfo<<"Writing Checkpoint \n"; + } + else + { + t += dt; + timestep++; + + gsField<> solution(mp,deformation); + if (timestep % plotmod==0 && plot) + { + // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + std::string fileName = "./output/solution" + util::to_string(timestep); + gsWriteParaview<>(solution, fileName, 500); + fileName = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileName,t,".vts"); + } + if (write) + { + solution.patch(0).eval_into(points,pointDataMatrix); + otherDataMatrix< +#include +#include +#include +// #include + +#include +#include + + +#ifdef gsStructuralAnalysis_ENABLED +#include +#include +#include +#include +#include +#include +#include +#endif + +#ifdef gsKLShell_ENABLED +#include +#include +#include +#include +#include +#endif + + +using namespace gismo; + +// template +// T ProjectForceMesh( const gsMatrix<> & forceControlPoints, +// const gsMultiBasis<> & forceBasis, +// const gsMultiBasis<> & geometryBasis, +// const index_t & targetDim, +// gsMatrix & results) +// { +// gsExprAssembler<> A(1,1); + +// gsMatrix solVector; + +// // Set the discretization space +// space u = A.getSpace(forceBasis, targetDim); + +// solution sol = A.getSolution(u, solVector); +// auto f = A.getCoordinates(u); + +// u.setup(0); +// A.initSystem(); + + + +// gsMatrix<> matrixForce, matrixGeometry; +// } + + +int main(int argc, char *argv[]) +{ + bool plot = false; + bool get_readTime = false; + bool get_writeTime = false; + index_t plotmod = 1; + index_t numRefine = 1; + index_t numElevate = 0; + int method = 3; // 1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson, 6 RK4 + + real_t rho = 3000; + real_t E = 4e6; + real_t nu = 0.3; + + //! [Parse command line] + std::string precice_config("../precice_config.xml"); + + gsCmdLine cmd("Coupled heat equation using 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("readTime", "Get the read time", get_readTime); + cmd.addSwitch("writeTime", "Get the write time", get_writeTime); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + GISMO_ASSERT(gsFileManager::fileExists(precice_config),"No precice config file has been defined"); + //! [Read input file] + std::string participantName = "Solid"; + gsPreCICE participant(participantName, precice_config); + + //Generate domain + gsMultiPatch<> patches; + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,0.0,0.5,1.0)); + + //Embed the 2D geometry to 3D + gsMultiPatch<> solutions; + patches.addAutoBoundaries(); + patches.embed(3); + + // Create bases + // p-refine + if (numElevate!=0) + patches.degreeElevate(numElevate); + + // h-refine + for (int r =0; r < numRefine; ++r) + patches.uniformRefine(); + + // Create bases + gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) + + gsInfo << "Patches: "<< patches.nPatches() <<", degree: "<< bases.minCwiseDegree() <<"\n"; + + + /* + * Data initialization + * + * This participant receives mesh information (knot vector, control points) from the solid, + * and it creates its own mesh based on the information. + * And writes pressure, reads displacements from the solid participant. + * The follow meshes and data are made available: + * + * - Meshes: + * + GeometryControlPointMesh: This mesh contains the control points as mesh vertices. + * + GeometryKnotMesh: This mesh contains the knots as mesh vertices. + * + ForceKnotMesh: This mesh contains the knots as mesh vertices for force mesh. + * + ForceControlPointMesh: This mesh contains the control points as mesh vertices for force mesh. + * + * + * - Data: + * + GeometryControlPointData: This data is defined on the GeometryControlPointMesh and stores the displacement of the control points + * + GeometryKnotData + * + ForceControlPointData: This data is defined on the ForceControlPointMesh and stores the change of pressure/stress + * + ForceKnotData + */ + + std::string GeometryKnotMesh = "GeometryKnotMesh"; + std::string GeometryControlPointMesh = "GeometryControlPointMesh"; + std::string ForceKnotMesh = "ForceKnotMesh"; + std::string ForceControlPointMesh = "ForceControlPointMesh"; + + // std::string StressData = "StressData"; + std::string GeometryControlPointData = "GeometryControlPointData"; + // std::string GeometryKnotData = "GeometryKnotData"; + // std::string ForceKnotData = "ForceKnotMeshData"; + std::string ForceControlPointData = "ForceControlPointData"; + + gsMatrix<> bbox(3,2); + bbox << -1e300, 1e300, // X dimension limits + -1e300, 1e300, // Y dimension limits + -1e300, 1e300; // Z dimension limits + bbox.transposeInPlace(); + bbox.resize(1,bbox.rows()*bbox.cols()); + participant.setMeshAccessRegion(ForceControlPointMesh, bbox); + + // Setup geometry control point mesh + gsVector geometryControlPointIDs; + gsMatrix<> geometryControlPoints = patches.patch(0).coefs().transpose(); + participant.addMesh(GeometryControlPointMesh, geometryControlPoints, geometryControlPointIDs); + + + // Setup the geometry knot mesh + gsMatrix<> geometryKnots = knotsToMatrix(bases.basis(0)); + + + participant.addMesh(GeometryKnotMesh,geometryKnots); + + real_t precice_dt = participant.initialize(); + + //G[et the force mesh from the coupling interface + gsVector forceKnotIDs; + gsMatrix<> forceKnots; + participant.getMeshVertexIDsAndCoordinates(ForceKnotMesh,forceKnotIDs,forceKnots); + + gsBasis<> * basis = knotMatrixToBasis(forceKnots).get(); + + gsVector forceControlPointIDs; + gsMatrix<> forceControlPoints; + participant.getMeshVertexIDsAndCoordinates(ForceControlPointMesh, forceControlPointIDs,forceControlPoints); + + + // // Step 2: Regenerate the geometry + gsMultiPatch<> forceMesh; //Geometry object belongs to gsFunctionSet + forceMesh.addPatch(give(basis->makeGeometry(forceControlPoints.transpose()))); + + + // Define boundary condition for solid mesh + gsBoundaryConditions<> bcInfo; + //Bottom Side + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, nullptr, -1); + bcInfo.addCondition(0, boundary::south, condition_type::clamped, nullptr, 2); + + // Assign geometry map + bcInfo.setGeoMap(patches); + + // Surface force equals to the force mesh + +//------------------------------------------ gsKLShell Module ------------------------------------------------ + //Setup the material properties + gsFunctionExpr<> E_modulus(std::to_string(E),3); + gsFunctionExpr<> PoissonRatio(std::to_string(nu),3); + gsFunctionExpr<> Density(std::to_string(rho),3); + + gsOptionList options; + + + //Define thickness + real_t thickness = 0.1; + + gsFunctionExpr<> t(std::to_string(thickness),3); + + std::vector*> parameters(2); + parameters[0] = &E_modulus; + parameters[1] = &PoissonRatio; + + gsMaterialMatrixBase* 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>(patches, t, parameters, Density, options); + + // REMOVE IT LATER + + // forceMesh.embed(3); + + gsThinShellAssembler<3, real_t, true> assembler(patches, bases, bcInfo, forceMesh, materialMatrix); + gsOptionList assemblerOptions = options.wrapIntoGroup("Assembler"); + + assembler.assemble(); + // forceMesh.patch(0).coefs().setZero(); + // assembler.assemble(); + // gsDebugVar(assembler.rhs()); + + assembler.setOptions(assemblerOptions); + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + + + // Compute the mass matrix (since it is constant over time) + assembler.assembleMass(); + gsSparseMatrix<> M = assembler.massMatrix(); + assembler.assemble(); + gsSparseMatrix<> K = assembler.matrix(); + + + // Define the solution collection for Paraview + std::string dirname = "./output"; + + gsFileManager::mkdir(dirname); + gsParaviewCollection collection(dirname + "/solution"); + + // Time step + real_t dt = precice_dt; + real_t t_read = 0; + real_t t_write = 0; + + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&solutions](gsMatrix const &x, gsSparseMatrix & m) + { + // to do: add time dependency of forcing + // For the shell + ThinShellAssemblerStatus status; + assembler.constructSolution(x, solutions); + status = assembler.assembleMatrix(solutions); + m = assembler.matrix(); + return true; + }; + + + // Function for the Residual + gsStructuralAnalysisOps::TResidual_t Residual = [&assembler,&solutions](gsMatrix const &x, real_t t, gsVector & result) + { + //Add assemble vector JL + ThinShellAssemblerStatus status; + assembler.constructSolution(x,solutions); + status = assembler.assembleVector(solutions); + 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; }; + + 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) + 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-3); + timeIntegrator->options().setSwitch("Verbose",true); + + + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + // RHS of the projection + gsMatrix<> solVector; + solVector.setZero(assembler.numDofs(),1); + + // Assemble the RHS + 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); + + + real_t time = 0; + // Plot initial solution + if (plot) + { + gsMultiPatch<> solution; + gsVector<> displacements = U; + solution = assembler.constructDisplacement(displacements); + solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + gsField<> solField(patches,solution); + std::string fileName = dirname + "/solution" + util::to_string(timestep); + gsWriteParaview<>(solField, fileName, 500); + fileName = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileName,time,".vts"); + } + + gsMatrix<> points(2,1); + points.col(0)<<0.5,1; + + gsStructuralAnalysisOutput writer("./output/pointData.csv",points); + writer.init({"x","y","z"},{"time"}); // point1 - x, point1 - y, point1 - z, time + + gsMatrix<> pointDataMatrix; + gsMatrix<> otherDataMatrix(1,1); + + // gsDebugVar(dt); + + // 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; + gsVector<> displacements = U; + solution = assembler.constructDisplacement(displacements); + // write heat fluxes to interface + geometryControlPoints = solution.patch(0).coefs().transpose(); + participant.writeData(GeometryControlPointMesh,GeometryControlPointData,geometryControlPointIDs,geometryControlPoints); + + if (get_writeTime) + t_write +=participant.writeTime(); + + + // 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 the time step + // advance variables + time += dt; + timestep++; + + gsField<> solField(patches,solution); + if (timestep % plotmod==0 && plot) + { + // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + std::string fileName = dirname + "/solution" + util::to_string(timestep); + gsWriteParaview<>(solField, fileName, 500); + fileName = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileName,time,".vts"); + } + solution.patch(0).eval_into(points,pointDataMatrix); + otherDataMatrix< + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/perpendicular-flap-IGA-Spine/precice-config.xml b/examples/perpendicular-flap-IGA-Spine/precice-config.xml new file mode 100644 index 0000000..491fae2 --- /dev/null +++ b/examples/perpendicular-flap-IGA-Spine/precice-config.xml @@ -0,0 +1,72 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/perpendicular-flap-IGA-Spine/solid-gismo-IGA.cpp b/examples/perpendicular-flap-IGA-Spine/solid-gismo-IGA.cpp new file mode 100644 index 0000000..1901261 --- /dev/null +++ b/examples/perpendicular-flap-IGA-Spine/solid-gismo-IGA.cpp @@ -0,0 +1,463 @@ +/** @file flow-over-heated-plate.cpp + + @brief Heat equation participant for the PreCICE example "flow over heated plate" + + 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 (TU Delft, 2023-...) + + Give knots/control points for the boundary curve +*/ + +#include +#include +#include +#include +// #include + +#include +#include + +#ifdef gsStructuralAnalysis_ENABLED +#include +#include +#include +#include +#include +#include + +#include +#endif + +using namespace gismo; + +template class basis_type> +inline void getKnots(const gsBasis & source, std::vector> & tensorKnots) +{ + if ( const basis_type * basis = dynamic_cast*>(&source) ) + for (index_t d=0; d!=DIM; d++) + tensorKnots[d] = basis->knots(d); +} + +int main(int argc, char *argv[]) +{ + //! [Parse command line] + bool plot = false; + index_t plotmod = 1; + index_t numRefine = 1; + index_t numElevate = 0; + std::string precice_config; + 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); + 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; + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(-0.05,0.0,0.05,1.0)); + + // p-refine + if (numElevate!=0) + patches.degreeElevate(numElevate); + + // h-refine + for (int r =0; r < numRefine; ++r) + patches.uniformRefine(); + + // patches.embed(3); + // Create bases + gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) + + gsInfo << "Patches: "<< patches.nPatches() <<", degree: "<< bases.minCwiseDegree() <<"\n"; + + real_t rho = 3000; + real_t E = 4e6; + real_t nu = 0.3; + + // 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 + * + * + */ + 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 KnotMesh = "KnotMesh"; + std::string ControlPointMesh= "ControlPointMesh"; + std::string ForceMesh = "ForceMesh"; + + std::string ControlPointData= "ControlPointData"; + std::string ForceData = "ForceData"; + + + + // Step 3: initialize the participant + // 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()); + gsPreCICEFunction g_C(&participant,ForceMesh,ForceData,patches,patches.geoDim(),false); + // Add all BCs + // Coupling interface + // bcInfo.addCondition(0, boundary::north, condition_type::neumann , &g_C); + bcInfo.addCondition(0, couplingInterfaces[0], condition_type::neumann , &g_C, -1, true); + bcInfo.addCondition(0, couplingInterfaces[1], condition_type::neumann , &g_C, -1, true); + bcInfo.addCondition(0, couplingInterfaces[2], condition_type::neumann , &g_C, -1, true); + // bcInfo.addCondition(0, boundary::west, condition_type::neumann , &g_C); + // 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); + + std::vector::uPtr> boundaries(couplingInterfaces.size()); + gsMultiPatch<> boundaryGeometry; + + // Derive the basis for boundaries + + // Finish this for both knot vectors and control points + gsMultiBasis<> boundaryBases; + + size_t vector_size; + size_t N = 0; + + // Pre-calculate total size N to resize derived_knots just once. + for(index_t i= 0; i < couplingInterfaces.size();++i) + { + boundaries[i] = patches.patch(0).boundary(couplingInterfaces[i].side()); // Add boundary coefficients to a vector + auto& basis_temp = boundaries[i]->basis(); // Assuming this returns a pointer to gsBasis + boundaryBases.addBasis(&basis_temp); + N += knotsToVector(boundaryBases.basis(i)).size(); + } + + gsVector<> derived_knots; + derived_knots.resize(N); + + // boundaryBases.basis(i).knots().asMatrix() + + // static_cast>(boundaryBases.basis(i)).knots().asMatrix() + + size_t offset = 0; // Tracks the current fill position in derived_knots. + for (index_t i = 0; i < couplingInterfaces.size(); ++i) + { + // No need to call knotsToVector twice for the same basis. + const gsVector<> temp_knots = knotsToVector(boundaryBases.basis(i)); + derived_knots.segment(offset, temp_knots.size()) = temp_knots; + offset += temp_knots.size(); // Increment offset by the size of the current knot vector. + } + gsDebugVar(derived_knots); + + // Control Points + // In one matrix + // gsMatrix<> derived_control_points; + + size_t totalControlPoints = 0; + size_t controlpoint_offset = 0; + for (index_t i = 0; i < couplingInterfaces.size(); ++i) + { + totalControlPoints += boundaries[i]->coefs().rows(); + } + gsMatrix<> derived_control_points(totalControlPoints, boundaries[0]->coefs().cols()); + derived_control_points.setZero(); + + for (index_t i = 0; i < couplingInterfaces.size(); ++i) + { + const gsMatrix<> temp_control_points = boundaries[i]->coefs(); + for (index_t j = 0; j < temp_control_points.rows(); ++j) + { + derived_control_points.row(controlpoint_offset + j) << temp_control_points.row(j); + } + controlpoint_offset += temp_control_points.rows(); + } + + gsDebugVar(derived_control_points); + + // Step 1: write the meshes to PreCICE + // Step 1a: KnotMesh + // get the knots in a matrix, using the utility function knotsToMatrix + gsMatrix<> knots = knotsToMatrix(bases.basis(0)); // + participant.addMesh(KnotMesh,derived_knots.transpose()); + + // Step 1b: ControlPointMesh + // get the control points, in the format where every column is a control point + gsVector controlPointIDs; // needed for writing + gsMatrix<> controlPoints = derived_control_points; + // gsDebugVar(controlPoints); + + participant.addMesh(ControlPointMesh,derived_control_points.transpose(),controlPointIDs); + + // Step 1c: ForceMesh + // Get the quadrature nodes on the coupling interface + gsOptionList quadOptions = gsAssembler<>::defaultOptions(); + + // Get the quadrature points + gsMatrix<> quadPoints = gsQuadrature::getAllNodes(bases.basis(0),quadOptions,couplingInterfaces); + gsDebugVar(quadPoints.dim()); + participant.addMesh(ForceMesh,quadPoints); + + + // Step 2 (not needed???) + // Needed for direct mesh coupling + gsMatrix<> bbox(2,2); + bbox.col(0).setConstant(-1e300); + bbox.col(1).setConstant(1e300); + bbox.transposeInPlace(); + participant.setMeshAccessRegion(ForceMesh,bbox); + + + real_t precice_dt = participant.initialize(); + + + + // for (index_t i = 0; i < couplingInterfaces.size(); ++i) { + // derived_control_points.push_back(boundaries[i]->coefs().transpose()); + // } + // std::cout << "Control Points: " << derived_control_points[0]<< std::endl; + + + // gsDebugVar(boundaryBases.basis(0)); + // derived_knots.push_back(knotsToVector(boundaryBases.basis(0))); + + // gsDebugVar(derived_knots); + + + // Maybe also make the fluid part as well (for testing). + +//---------------------------------------------------------------------------------------------- + // // 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::hooke); + 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<> 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); + + // 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; + }; + + 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; }; + + 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) + 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-3); + timeIntegrator->options().setSwitch("Verbose",true); + + real_t time = 0; + // Plot initial solution + if (plot) + { + gsMultiPatch<> solution; + assembler.constructSolution(solVector,fixedDofs,solution); + + // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + 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(2,1); + points.col(0)<<0.5,1; + + gsStructuralAnalysisOutput writer("pointData.csv",points); + writer.init({"x","y"},{"time"}); // point1 - x, point1 - y, time + + gsMatrix<> pointDataMatrix; + gsMatrix<> otherDataMatrix(1,1); + + gsDebugVar(dt); + + // 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 + // solution.embed(3); + controlPoints = solution.patch(0).coefs().transpose(); + + participant.writeData(ControlPointMesh,ControlPointData,controlPointIDs,controlPoints); + + // 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 && 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"); + } + + solution.patch(0).eval_into(points,pointDataMatrix); + otherDataMatrix< +#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 = 1; + index_t loadCase = 0; + bool get_readTime = false; + bool get_writeTime = false; + std::string precice_config; + + gsCmdLine cmd("Flow over heated plate for PreCICE."); + 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("l","loadCase", "Load case: 0=constant load, 1='spring' load", loadCase); + cmd.addSwitch("readTime", "Get the read time", get_readTime); + cmd.addSwitch("writeTime", "Get the write time", get_writeTime); + 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; + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,0.0,0.5,1.0)); + + patches.addAutoBoundaries(); + patches.embed(3); + + // Generate discrete fluid mesh + gsMatrix<> bbox = patches.patch(0).support(); + + int numPoints = 100; + gsMatrix<> pointGrid = gsPointGrid<>(bbox, numPoints); + gsMatrix<> mesh = patches.patch(0).eval(pointGrid); + + // Embed the 2D geometry to 3D + + /* + * Initialize the preCICE participant + * + * + */ + std::string participantName = "Fluid"; + gsPreCICE participant(participantName, precice_config); + + /* + * Data initialization + * + * This participant creates its own mesh, and it writes and reads displacements and stresses on its mesh. + * The follow meshes and data are made available: + * + * - Meshes: + * + FluidMesh: This mesh contains the integration points as mesh vertices + * + * - Data: + * + DisplacementData: This data is defined on the FluidMesh and stores the displacement at the integration points + * + StressData: This data is defined on the FluidMesh and stores pressure/forces at the integration points + */ + std::string FluidMesh = "FluidMesh"; + std::string StressData = "StressData"; + std::string DisplacementData = "DisplacementData"; + + + // Step 1: FluidMesh + // Get the quadrature points + gsVector FluidMeshIDs; + gsDebugVar(mesh.dim()); + participant.addMesh(FluidMesh,mesh,FluidMeshIDs); + gsDebugVar(mesh); + + + // Step 2: initialize the participant + real_t precice_dt = participant.initialize(); + gsDebugVar("Initialize Fluid"); + + /* + * Collect the geometry + * + */ + + real_t t = 0, dt = precice_dt; + real_t t_read = 0; + real_t t_write = 0; + index_t timestep = 0; + // Define the solution collection for Paraview + gsParaviewCollection collection("./output/solution"); + + gsMatrix<> StressPointData(3,mesh.cols()); + + // Time integration loop + while (participant.isCouplingOngoing()) + { + if (participant.requiresWritingCheckpoint()) + gsInfo<<"Writing Checkpoint\n"; + + // Read control point displacements + gsMatrix<> meshPointDisplacements; + participant.readData(FluidMesh,DisplacementData,FluidMeshIDs,meshPointDisplacements); + + if (get_readTime) + t_read += participant.readTime(); + + + if (loadCase==0) + { + // Write data at the quadrature points + StressPointData.setZero(); + StressPointData.row(2).setConstant(-5e3); + } + else if (loadCase==1) + { + StressPointData.setZero(); + for (index_t k = 0; k!=meshPointDisplacements.cols(); k++) + { + StressPointData(2,k) = -meshPointDisplacements(2,k); + } + // impact loading at t==0 + if (t==0) + StressPointData.row(2).array() += -5e3; + } + else + GISMO_ERROR("Load case "< solution(mp,deformation); + if (timestep % plotmod==0 && plot) + { + // // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + // std::string fileName = "./output/solution" + util::to_string(timestep); + // gsWriteParaview<>(solution, fileName, 500); + + // fileName = "solution" + util::to_string(timestep) + "0"; + // collection.addTimestep(fileName,t,".vts"); + } + + // solution.patch(0).eval_into(points,pointDataMatrix); + // otherDataMatrix< + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/perpendicular-flap-vertex-vertex/solid-gismo/perpendicular-flap-solid-gismo.cpp b/examples/perpendicular-flap-vertex-vertex/solid-gismo/perpendicular-flap-solid-gismo.cpp new file mode 100644 index 0000000..6d12c4b --- /dev/null +++ b/examples/perpendicular-flap-vertex-vertex/solid-gismo/perpendicular-flap-solid-gismo.cpp @@ -0,0 +1,564 @@ +/** @file flow-over-heated-plate.cpp + + @brief Heat equation participant for the PreCICE example "flow over heated plate" + + 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/. + + Solid solver: quadrature point + Fluid solver: displacements on its control points + + Author(s): J.Li (TU Delft, 2023-...) +*/ + +#include +#include +#include +#include +#include +// #include + +#include +#include + + +#ifdef gsKLShell_ENABLED +#include +#include +#include +#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; + bool write = false; + bool get_readTime = false; + bool get_writeTime = false; + index_t plotmod = 1; + index_t numRefine = 2; + index_t numElevate = 1; + std::string precice_config; + int method = 3; // 1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson, 6 RK4 + + 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","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",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); + 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; + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(-0.05,0.0,0.05,1.0)); + + // Embed the 2D geometry to 3D + gsMultiPatch<> solutions; + patches.addAutoBoundaries(); + patches.embed(3); + // patches.embed(3); + // source function, rhs + gsConstantFunction<> g(0.,0.,3); + + // source function, rhs + gsConstantFunction<> gravity(0.,0.,3); + + // Create bases + // p-refine + if (numElevate!=0) + patches.degreeElevate(numElevate); + + // h-refine + for (int r =0; r < numRefine; ++r) + patches.uniformRefine(); + + // Create bases + gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) + + gsInfo << "Patches: "<< patches.nPatches() <<", degree: "<< bases.minCwiseDegree() <<"\n"; + + real_t rho = 3000; + real_t E = 4e6; + real_t nu = 0.3; + // real_t nu = 0.3; + real_t mu = E / (2.0 * (1.0 + nu)); + real_t lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); + + /* + * Initialize the preCICE participant + * + * + */ + std::string participantName = "Solid"; + gsPreCICE participant(participantName, precice_config); + + // ---------------------------------------------------------------------------------------------- + // typedef gsExprAssembler<>::geometryMap geometryMap; + // typedef gsExprAssembler<>::space space; + // typedef gsExprAssembler<>::solution solution; + + // gsExprAssembler<> A_quad(1,1); + + // gsInfo<<"Active options:\n"<< A.options() <<"\n"; + + // A.setIntegrationElements(bases); + + // gsExprEvaluator<> ev(A); + + + + + + /* + * Data initialization + * + * This participant creates its own mesh, and it writes and reads displacements and stresses on its mesh. + * The follow meshes and data are made available: + * + * - Meshes: + * + SolidMesh: This mesh contains the integration points as mesh vertices + * + * - Data: + * + DisplacementData: This data is defined on the SolidMesh and stores the displacement at the integration points + * + StressData: This data is defined on the SolidMesh and stores pressure/forces at the integration points + */ + std::string SolidMesh = "Solid-Mesh"; + std::string StressData = "Force"; + std::string DisplacementData = "Displacement"; + // std::string ForceMesh = "Fluid-Mesh"; + + std::vector couplingInterfaces(3); + couplingInterfaces[0] = patchSide(0,boundary::west); + couplingInterfaces[1] = patchSide(0,boundary::north); + couplingInterfaces[2] = patchSide(0,boundary::east); + + // gsOptionList quadOptions = A_quad.options(); + + + // Step 1: SolidMesh + // Get the quadrature nodes on the coupling interface + gsOptionList quadOptions = gsExprAssembler<>::defaultOptions(); + + // Get the quadrature points + gsVector quadPointIDs; + // gsMatrix<> quadPoints = gsQuadrature::getAllNodes(bases.basis(0),quadOptions,couplingInterfaces); + // gsMatrix<> datapoints = patches.patch(0).eval(quadPoints); + // gsMatrix<> quadPointsUpdate(2,13); + + // Set the first point to (0, 0) + // quadPointsUpdate(0, 0) = 0; + // quadPointsUpdate(1, 0) = 0; + + // index_t quadSize = 0; + + // for (std::vector::const_iterator it = couplingInterfaces.begin(); it!=couplingInterfaces.end(); it++) + // { + // // Get a domain iterator on the coupling interface + // typename gsBasis::domainIter domIt = bases.basis(it->patch).makeDomainIterator(it->side()); + + // // First obtain the size of all quadrature points + // typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT + // for (; domIt->good(); domIt->next() ) + // { + // QuRule = gsQuadrature::getPtr(bases.basis(it->patch), quadOptions,it->side().direction()); + // quadSize+=QuRule->numNodes(); + // } + // } + // quadPointsUpdate.block(0, 1, 2, 12) = datapoints.block(0, 0, 2, 12); + + // gsDebugVar(quadPointsUpdate); + + + + + // Initialize parametric coordinates + gsMatrix<> quadPoints = gsQuadrature::getAllNodes(bases.basis(0),quadOptions); + + + // gsMatrix<> quadPoints(patches.domainDim(),quadSize); + // quadPoints.setZero(); + // // Initialize physical coordinates + // gsMatrix<> datapoints(patches.targetDim(),quadSize); + // datapoints.setZero(); + + // Grab all quadrature points + index_t offset = 0; + // Set the dimension of the points + gsMatrix<> nodes; + // Start iteration over elements + gsVector<> tmp; + + + // for (std::vector::const_iterator it = couplingInterfaces.begin(); it!=couplingInterfaces.end(); it++) + // { + // // Get a domain iterator on the coupling interface + // typename gsBasis::domainIter domIt = bases.basis(it->patch).makeDomainIterator(it->side()); + // typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT + // for (domIt->reset(); domIt->good(); domIt->next()) + // { + // QuRule = gsQuadrature::getPtr(bases.basis(it->patch), quadOptions,it->side().direction()); + // // Map the Quadrature rule to the element + // QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(), + // nodes, tmp); + // quadPoints.block(0,offset,patches.domainDim(),QuRule->numNodes()) = nodes; + + // gsMatrix<> tmp2; + // patches.patch(it->patch).eval_into(nodes,tmp2); + // datapoints.block(0,offset,patches.targetDim(),QuRule->numNodes()) = patches.patch(it->patch).eval(nodes); + // offset += QuRule->numNodes(); + // } + // } + + + gsDebugVar(quadPoints); + + + + participant.addMesh(SolidMesh,quadPoints,quadPointIDs); //Set the vertices to be datapoints (quadpoints in physical domain) + + + + gsMatrix<> quadPointsData(3, quadPoints.cols()); + + + // gsDebugVar(quadPoints); + // gsDebugVar(quadPoints.dim()); + quadPointsData.setZero(); + + + + // Step 2: initialize the participant + real_t precice_dt = participant.initialize(); + + gsMultiPatch<> forceMesh; + + gsMultiBasis<> forceBases; + std::vector::uPtr> couplingBoundaries(couplingInterfaces.size()); + couplingBoundaries[0] = patches.patch(0).boundary(couplingInterfaces[0].side()); + couplingBoundaries[1] = patches.patch(0).boundary(couplingInterfaces[1].side()); + couplingBoundaries[2] = patches.patch(0).boundary(couplingInterfaces[2].side()); + + + //--------------------------Option: using spline to represent the stress field------------------- + forceMesh.addPatch(patches.patch(0).boundary(couplingInterfaces[0].side())); + forceMesh.addPatch(patches.patch(0).boundary(couplingInterfaces[1].side())); + forceMesh.addPatch(patches.patch(0).boundary(couplingInterfaces[2].side())); + + gsDebugVar(forceMesh.basis(0).size()); + gsDebugVar(forceMesh.basis(1).size()); + + index_t N_cps = forceMesh.basis(0).size(); //Number of control points in y direction + index_t M_cps = forceMesh.basis(1).size(); //Number of control points in y direction + + + + + // gsDebugVar(forceMesh.basis(0).eval(quadPoints)); + + + //Collect the geometry + + // participant.getMeshVertexIDsAndCoordinates(ForceMesh,quadPointIDs,quadPoints); + + + + +//---------------------------------------------------------------------------------------------- + + // Define boundary conditions + gsBoundaryConditions<> bcInfo; + + // Bottom side + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, nullptr, -1); + bcInfo.addCondition(0, boundary::south, condition_type::clamped, nullptr, 2); + + // // Top side + // bcInfo.addCondition(0, boundary::west, condition_type::neumann, &forceMesh.patch(0)); + // bcInfo.addCondition(0, boundary::north, condition_type::neumann, &forceMesh.patch(1)); + // bcInfo.addCondition(0, boundary::east, condition_type::neumann, &forceMesh.patch(2)); + + // Assign geometry map + bcInfo.setGeoMap(patches); + gsLookupFunction surfForce(quadPoints, quadPointsData); + + gsDebugVar(quadPoints); + gsDebugVar(quadPointsData); + + // Set up the material matrices + gsFunctionExpr<> E_modulus(std::to_string(E),3); + gsFunctionExpr<> PoissonRatio(std::to_string(nu),3); + gsFunctionExpr<> Density(std::to_string(rho),3); + + //Define thickness + real_t thickness = 0.1; + + gsFunctionExpr<> t(std::to_string(thickness),3); + + std::vector*> parameters(2); + parameters[0] = &E_modulus; + parameters[1] = &PoissonRatio; + + gsOptionList options; + + gsMaterialMatrixBase* 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>(patches, t, parameters, Density, options); + + + gsThinShellAssembler<3, real_t, false> assembler(patches, bases, bcInfo, surfForce, materialMatrix); + gsOptionList assemblerOptions = quadOptions.wrapIntoGroup("ExprAssembler"); + assembler.setOptions(assemblerOptions); + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + + // Compute the mass matrix (since it is constant over time) + assembler.assembleMass(); + gsSparseMatrix<> M = assembler.massMatrix(); + gsInfo << "Got here \n"; + assembler.assemble(); + gsSparseMatrix<> K = assembler.matrix(); + + + // Define the solution collection for Paraview + gsFileManager::mkdir(dirname); + gsParaviewCollection collection(dirname + "/solution"); + + // Time step + real_t dt = precice_dt; + + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&solutions](gsMatrix const &x, gsSparseMatrix & m) + { + // to do: add time dependency of forcing + // For the shell + ThinShellAssemblerStatus status; + assembler.constructSolution(x, solutions); + status = assembler.assembleMatrix(solutions); + m = assembler.matrix(); + return true; + }; + + + // Function for the Residual + gsStructuralAnalysisOps::TResidual_t Residual = [&assembler,&solutions](gsMatrix const &x, real_t t, gsVector & result) + { + //Add assemble vector JL + ThinShellAssemblerStatus status; + assembler.constructSolution(x,solutions); + status = assembler.assembleVector(solutions); + 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; }; + + 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) + 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-3); + timeIntegrator->options().setSwitch("Verbose",true); + + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + // RHS of the projection + gsMatrix<> solVector; + solVector.setZero(assembler.numDofs(),1); + + // Assemble the RHS + 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); + + + real_t time = 0; + real_t t_read = 0; + real_t t_write = 0; + + // // Plot initial solution + // if (plot) + // { + // gsMultiPatch<> solution; + // gsVector<> displacements = U; + // solution = assembler.constructDisplacement(displacements); + + // // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + // gsField<> solField(patches,solution); + + // std::string fileName = dirname + "/solution" + util::to_string(timestep); + // gsWriteParaview<>(solField, fileName, 500); + // fileName = "solution" + util::to_string(timestep) + "0"; + // collection.addTimestep(fileName,time,".vts"); + // } + + gsMatrix<> points(2,1); + points.col(0)<<0.5,1; + + + gsStructuralAnalysisOutput writer("./output/pointData.csv",points); + writer.init({"x","y","z"},{"time"}); // point1 - x, point1 - y, point1 - z, time + + gsMatrix<> pointDataMatrix; + gsMatrix<> otherDataMatrix(1,1); + + gsDebugVar("Got here"); + + // 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|| = "< tempData(2,quadPointsData.cols()); + dt = std::min(dt,precice_dt); + participant.readData(SolidMesh,StressData, quadPointIDs, tempData); + + quadPointsData.row(0) = tempData.row(0).array() ; + quadPointsData.row(1) = tempData.row(1).array() ; + + gsDebugVar(quadPointsData); + assembler.assemble(); + F = assembler.rhs(); + + // 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 + + gsMultiPatch<> solution; + gsVector<> displacements = U; + solution = assembler.constructDisplacement(displacements); + // Information to fluid needs to be 2D + gsMatrix<> tempDisplacement(2, quadPoints.cols()); + gsMatrix<> quadPointsDisplacements = solution.patch(0).eval(quadPoints); + tempDisplacement.row(0) = quadPointsDisplacements.row(0); + tempDisplacement.row(1) = quadPointsDisplacements.row(1); + gsDebugVar(quadPointsDisplacements); + participant.writeData(SolidMesh,DisplacementData,quadPointIDs,tempDisplacement); + + if (get_writeTime) + t_write +=participant.writeTime(); + + // 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 && plot) + { + // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + std::string fileName = dirname + "/solution" + util::to_string(timestep); + gsWriteParaview<>(solField, fileName, 500); + fileName = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileName,time,".vts"); + } + solution.patch(0).eval_into(points,pointDataMatrix); + if (get_readTime) + t_read += participant.readTime(); + otherDataMatrix< +#include +#include +//#include +#include +// #include + +#include +#include + + +#ifdef gsKLShell_ENABLED +#include +#include +#include +#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; + bool write = false; + bool get_readTime = false; + bool get_writeTime = false; + index_t plotmod = 1; + index_t numRefine = 1; + index_t numElevate = 0; + std::string precice_config; + int method = 3; // 1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson, 6 RK4 + + 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","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",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); + 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; + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,0.0,0.5,1.0)); + + // Embed the 2D geometry to 3D + gsMultiPatch<> solutions; + patches.addAutoBoundaries(); + patches.embed(3); + + + + + // Create bases + // p-refine + if (numElevate!=0) + patches.degreeElevate(numElevate); + + // h-refine + for (int r =0; r < numRefine; ++r) + patches.uniformRefine(); + + // Create bases + gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) + + gsInfo << "Patches: "<< patches.nPatches() <<", degree: "<< bases.minCwiseDegree() <<"\n"; + + real_t rho = 3000; + real_t E = 4e6; + real_t nu = 0.3; + // real_t nu = 0.3; + // real_t mu = E / (2.0 * (1.0 + nu)); + // real_t lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); + + /* + * Initialize the preCICE participant + * + * + */ + std::string participantName = "Solid"; + gsPreCICE participant(participantName, precice_config); + + + + + /* + * Data initialization + * + * This participant creates its own mesh, and it writes and reads displacements and stresses on its mesh. + * The follow meshes and data are made available: + * + * - Meshes: + * + SolidMesh: This mesh contains the integration points as mesh vertices + * + * - Data: + * + DisplacementData: This data is defined on the SolidMesh and stores the displacement at the integration points + * + StressData: This data is defined on the SolidMesh and stores pressure/forces at the integration points + */ + std::string SolidMesh = "SolidMesh"; + std::string StressData = "StressData"; + std::string DisplacementData = "DisplacementData"; + + + // Step 1: SolidMesh + // Get the quadrature nodes on the coupling interface + gsOptionList quadOptions = gsExprAssembler<>::defaultOptions(); + + // Get the quadrature points + gsVector quadPointIDs; + gsMatrix<> quadPoints = gsQuadrature::getAllNodes(bases.basis(0),quadOptions); + participant.addMesh(SolidMesh,quadPoints,quadPointIDs); + gsMatrix<> quadPointsData(3, quadPoints.cols()); + quadPointsData.setZero(); + + + // Step 2: initialize the participant + real_t precice_dt = participant.initialize(); + +//---------------------------------------------------------------------------------------------- + + // Define boundary conditions + gsBoundaryConditions<> bcInfo; + + // Bottom side + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, nullptr, -1); + bcInfo.addCondition(0, boundary::south, condition_type::clamped, nullptr, 2); + + // Assign geometry map + bcInfo.setGeoMap(patches); + + // Surface force + // gsPreCICEFunction surfForce(&participant,SolidMesh,StressData,patches,patches.parDim(),patches.geoDim(),false); + gsLookupFunction surfForce(quadPoints, quadPointsData); + + + + +//------------------------------------------IF gsElasticity is Used------------------------------ + // // 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::hooke); + // assembler.assemble(); + + + // gsMatrix Minv; + // gsSparseMatrix<> M = massAssembler.matrix(); + // gsSparseMatrix<> K = assembler.matrix(); + // gsSparseMatrix<> K_T; +//------------------------------------------IF gsElasticity is Used------------------------------ + + + // Set up the material matrices + gsFunctionExpr<> E_modulus(std::to_string(E),3); + gsFunctionExpr<> PoissonRatio(std::to_string(nu),3); + gsFunctionExpr<> Density(std::to_string(rho),3); + + //Define thickness + real_t thickness = 0.1; + + gsFunctionExpr<> t(std::to_string(thickness),3); + + std::vector*> parameters(2); + parameters[0] = &E_modulus; + parameters[1] = &PoissonRatio; + + gsOptionList options; + + gsMaterialMatrixBase* 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>(patches, t, parameters, Density, options); + + gsThinShellAssembler<3, real_t, true> assembler(patches, bases, bcInfo, surfForce, materialMatrix); + gsOptionList assemblerOptions = quadOptions.wrapIntoGroup("ExprAssembler"); + assembler.setOptions(assemblerOptions); + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + + // Compute the mass matrix (since it is constant over time) + assembler.assembleMass(); + gsSparseMatrix<> M = assembler.massMatrix(); + assembler.assemble(); + gsSparseMatrix<> K = assembler.matrix(); + + + // Define the solution collection for Paraview + gsFileManager::mkdir(dirname); + gsParaviewCollection collection(dirname + "/solution"); + + // Time step + real_t dt = precice_dt; + + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&solutions](gsMatrix const &x, gsSparseMatrix & m) + { + // to do: add time dependency of forcing + // For the shell + ThinShellAssemblerStatus status; + assembler.constructSolution(x, solutions); + status = assembler.assembleMatrix(solutions); + m = assembler.matrix(); + return true; + }; + + + // Function for the Residual + gsStructuralAnalysisOps::TResidual_t Residual = [&assembler,&solutions](gsMatrix const &x, real_t t, gsVector & result) + { + //Add assemble vector JL + ThinShellAssemblerStatus status; + assembler.constructSolution(x,solutions); + status = assembler.assembleVector(solutions); + 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; }; + + 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) + 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-3); + timeIntegrator->options().setSwitch("Verbose",true); + + + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + // RHS of the projection + gsMatrix<> solVector; + solVector.setZero(assembler.numDofs(),1); + + // Assemble the RHS + 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); + + + real_t time = 0; + real_t t_read = 0; + real_t t_write = 0; + // Plot initial solution + if (plot) + { + gsMultiPatch<> solution; + gsVector<> displacements = U; + solution = assembler.constructDisplacement(displacements); + + // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + gsField<> solField(patches,solution); + + std::string fileName = dirname + "/solution" + util::to_string(timestep); + gsWriteParaview<>(solField, fileName, 500); + fileName = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileName,time,".vts"); + } + + gsMatrix<> points(2,1); + points.col(0)<<0.5,1; + + + gsStructuralAnalysisOutput writer("./output/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; + gsVector<> displacements = U; + solution = assembler.constructDisplacement(displacements); + gsMatrix<> quadPointsDisplacements = solution.patch(0).eval(quadPoints); + participant.writeData(SolidMesh,DisplacementData,quadPointIDs,quadPointsDisplacements); + + if (get_writeTime) + t_write +=participant.writeTime(); + + // 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 && plot) + { + // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + std::string fileName = dirname + "/solution" + util::to_string(timestep); + gsWriteParaview<>(solField, fileName, 500); + fileName = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileName,time,".vts"); + } + solution.patch(0).eval_into(points,pointDataMatrix); + if (get_readTime) + t_read += participant.readTime(); + otherDataMatrix< class gsLookupFunction : public gsFunction diff --git a/src/gsPreCICE.h b/src/gsPreCICE.h index af48abd..127c0a5 100644 --- a/src/gsPreCICE.h +++ b/src/gsPreCICE.h @@ -8,9 +8,9 @@ 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): H.M. Verhelst (TU Delft, 2019-...), J.Li (TU Delft, 2023-...) + Author(s): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) - Summary of Functions: + @details Summary of Functions: 1. Constructors: - gsPreCICE() DEFAULT - gsPreCICE(std::string participantName, std::string configurationFileName): diff --git a/src/gsPreCICEFunction.h b/src/gsPreCICEFunction.h index a9613ae..bdfcb94 100644 --- a/src/gsPreCICEFunction.h +++ b/src/gsPreCICEFunction.h @@ -8,7 +8,7 @@ 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): H.M. Verhelst (TU Delft, 2019-...) + Author(s): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) */ #pragma once @@ -22,7 +22,7 @@ namespace gismo /** * @brief Class defining a gsFunction that reads from the precice::SolverInterface - * + * @usage Combine with gsElasticity module to update boundary conditions. * @tparam T Number format */ template diff --git a/src/gsPreCICEUtils.h b/src/gsPreCICEUtils.h index 6116897..391265f 100644 --- a/src/gsPreCICEUtils.h +++ b/src/gsPreCICEUtils.h @@ -8,7 +8,7 @@ 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): H.M. Verhelst (TU Delft, 2019-...) + Author(s): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) */ #pragma once