From e433eb57d2560cbe9c05c7b9a3b0712bc729e8cf Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Thu, 12 Sep 2024 10:57:41 +0200 Subject: [PATCH 01/29] add torsion computation for annulus and cylinder --- filedata/pde/annulus_4p.xml | 6 +- filedata/pde/annulus_4p_TFT.xml | 6 +- filedata/pde/annulus_4p_hyperelastic.xml | 6 +- filedata/pde/cylinder_4p.xml | 6 +- filedata/pde/cylinder_4p_TFT.xml | 6 +- .../static_shell_multipatch_XML_torsion.cpp | 717 ++++++++++++++++++ 6 files changed, 742 insertions(+), 5 deletions(-) create mode 100644 solvers/static_shell_multipatch_XML_torsion.cpp diff --git a/filedata/pde/annulus_4p.xml b/filedata/pde/annulus_4p.xml index 8beb4dc..844c7cd 100644 --- a/filedata/pde/annulus_4p.xml +++ b/filedata/pde/annulus_4p.xml @@ -34,7 +34,11 @@ - 0 + + 0 + 0 + 0 + sqrt(x^2+y^2)*(cos(atan2(y,x) + pi/2)- cos(atan2(y,x))) sqrt(x^2+y^2)*(sin(atan2(y,x) + pi/2)- sin(atan2(y,x))) 125e-3 diff --git a/filedata/pde/annulus_4p_TFT.xml b/filedata/pde/annulus_4p_TFT.xml index 8ad7d3d..d0ddec6 100644 --- a/filedata/pde/annulus_4p_TFT.xml +++ b/filedata/pde/annulus_4p_TFT.xml @@ -34,7 +34,11 @@ - 0 + + 0 + 0 + 0 + sqrt(x^2+y^2)*(cos(atan2(y,x) + pi/2)- cos(atan2(y,x))) sqrt(x^2+y^2)*(sin(atan2(y,x) + pi/2)- sin(atan2(y,x))) 125e-3 diff --git a/filedata/pde/annulus_4p_hyperelastic.xml b/filedata/pde/annulus_4p_hyperelastic.xml index 9b691c4..c03a724 100644 --- a/filedata/pde/annulus_4p_hyperelastic.xml +++ b/filedata/pde/annulus_4p_hyperelastic.xml @@ -34,7 +34,11 @@ - 0 + + 0 + 0 + 0 + sqrt(x^2+y^2)*(cos(atan2(y,x) + pi/2)- cos(atan2(y,x))) sqrt(x^2+y^2)*(sin(atan2(y,x) + pi/2)- sin(atan2(y,x))) 125e-3 diff --git a/filedata/pde/cylinder_4p.xml b/filedata/pde/cylinder_4p.xml index 807eb80..d08f84c 100644 --- a/filedata/pde/cylinder_4p.xml +++ b/filedata/pde/cylinder_4p.xml @@ -34,7 +34,11 @@ - 0 + + 0 + 0 + 0 + sqrt(x^2+y^2)*(cos(atan2(y,x) + pi/2)- cos(atan2(y,x))) sqrt(x^2+y^2)*(sin(atan2(y,x) + pi/2)- sin(atan2(y,x))) 1 diff --git a/filedata/pde/cylinder_4p_TFT.xml b/filedata/pde/cylinder_4p_TFT.xml index 5703b25..a5d1553 100644 --- a/filedata/pde/cylinder_4p_TFT.xml +++ b/filedata/pde/cylinder_4p_TFT.xml @@ -34,7 +34,11 @@ - 0 + + 0 + 0 + 0 + sqrt(x^2+y^2)*(cos(atan2(y,x) + pi/2)- cos(atan2(y,x))) sqrt(x^2+y^2)*(sin(atan2(y,x) + pi/2)- sin(atan2(y,x))) 1 diff --git a/solvers/static_shell_multipatch_XML_torsion.cpp b/solvers/static_shell_multipatch_XML_torsion.cpp new file mode 100644 index 0000000..c639613 --- /dev/null +++ b/solvers/static_shell_multipatch_XML_torsion.cpp @@ -0,0 +1,717 @@ +/** @file static_shell_multipatch_XML.cpp + + @brief Blackbox solver for static shell analysis using unstructured splines + + 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 (2019-..., TU Delft) +*/ + +#include + +#ifdef gsKLShell_ENABLED +#include +#include +#include +#endif + +#include +#include +#include + +#include + +#ifdef gsUnstructuredSplines_ENABLED +#include +#include +#include +#endif + +#include + +using namespace gismo; + +#ifdef gsKLShell_ENABLED +#ifdef gsUnstructuredSplines_ENABLED +int main (int argc, char** argv) +{ + // Input options + int numElevate = 0; + int numRefine = 1; + bool plot = false; + bool write = false; + bool stress = false; + bool DR = false; + bool NR = false; + + bool membrane = false; + + index_t method = 0; + + real_t perturb = 0; + + // Arc length method options + + std::string bvp; + + std::string dirname = "."; + + std::string wn("data.csv"); + + gsCmdLine cmd("Shell static solver for multipatches."); + + cmd.addInt("r","hRefine", "Number of dyadic h-refinement (bisection) steps to perform before solving", numRefine); + cmd.addInt("e","degreeElevation", "Number of degree elevation steps to perform on the Geometry's basis before solving", numElevate); + cmd.addInt("m","method", "Smoothing method to use: 0: smoothInterfaces, 1: Almost C1, 2: D-Patch", method); + + cmd.addReal( "P" , "perturb" , "Perturb the z coordinate of the refined geometry", perturb ); + + cmd.addSwitch("membrane", "Membrane element", membrane); + + cmd.addSwitch("plot", "Plot result in ParaView format", plot); + cmd.addSwitch("write", "Write data to file", write); + cmd.addSwitch("stress", "Plot stress in ParaView format", stress); + cmd.addSwitch("DR", "Use Dynamic Relaxation", DR); + cmd.addSwitch("NR", "Use Newton Raphson", NR); + + cmd.addString("i","inputFile", "Input file", bvp); + cmd.addString("o","outputDir", "Output directory", dirname); + + + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + // ![Initialize members] + gsMultiPatch<> mp; + // Boundary conditions + gsBoundaryConditions<> BCs; + + gsFunctionExpr<> forceFun; + gsFunctionExpr<> pressFun; + + gsPointLoads pLoads = gsPointLoads(); + gsMatrix<> points,loads; + gsMatrix pid_ploads; + + gsOptionList DROptions, NROptions, assemblerOptions; + // ![Initialize members] + + // ![Read data] + gsFileData fd(bvp); + if (fd.hasAny>()) + { + gsInfo<<"Reading geometry from "<(path + geomFileName,mp); + } + gsInfo<<"Finished\n"; + + // p-refine + for (size_t p=0; p!=mp.nPatches(); p++) + { + for(index_t i = 0; i< numElevate; ++i) + { + if (dynamic_cast * >(&mp.patch(p))) + { + gsWarn<<"Degree elevation applied"<<"\n"; + mp.patch(p).degreeElevate(); // Elevate the degree + } + else + mp.patch(p).degreeIncrease(); // Elevate the degree + } + + // h-refine + for(index_t i = 0; i< numRefine; ++i) + mp.patch(p).uniformRefine(); + } + + if (perturb!=0) + for (size_t p=0; p!=mp.nPatches(); p++) + { + mp.patch(p).coefs().col(2).setRandom(); + mp.patch(p).coefs().col(2) *= perturb; + } + + gsMultiBasis<> dbasis(mp,true); + gsInfo<<"Basis (patch 0): "<< mp.patch(0).basis() << "\n"; + if (gsTensorNurbsBasis<2,real_t> * basis = dynamic_cast * >(&mp.basis(0))) + for (index_t dim = 0; dim!=2; dim++) + gsDebug<<"dir "<knots(dim).asMatrix()<<"\n"; + + gsInfo<<"Looking for material matrices ...\n"; + gsMaterialMatrixContainer materialMatrixContainer; + gsMaterialMatrixBase * materialMatrix; + if (fd.hasAny>()) + { + gsInfo<<"Reading material matrix container (ID=11) ..."; + fd.getId(11,materialMatrixContainer); + } + else + { + gsInfo<<"Reading material matrix (ID=10) ..."; + materialMatrix = fd.getId>(10).release(); + for (size_t p = 0; p!=mp.nPatches(); p++) + materialMatrixContainer.add(materialMatrix); + } + gsInfo<<"Finished\n"; + gsInfo<<"Finished\n"; + + gsInfo<<"Reading boundary conditions (ID=20) ..."; + fd.getId(20,BCs); + gsInfo<<"Finished\n"; + + gsInfo<<"Reading force function (ID=21) ..."; + fd.getId(21,forceFun); + gsInfo<<"Finished\n"; + bool pressure = false; + if ( fd.hasId(22) ) + { + gsInfo<<"Reading pressure function (ID=22) ..."; + pressure = true; + fd.getId(22,pressFun); + gsInfo<<"Finished\n"; + } + // Point loads + gsInfo<<"Reading point load locations (ID=30) ..."; + if ( fd.hasId(30) ) fd.getId(30,points); + gsInfo<<"Finished\n"; + gsInfo<<"Reading point load vectors (ID=31) ..."; + if ( fd.hasId(31) ) fd.getId(31,loads); + gsInfo<<"Finished\n"; + gsInfo<<"Reading point load patch indices (ID=32) ..."; + if ( fd.hasId(32) ) fd.getId(32,pid_ploads); + gsInfo<<"Finished\n"; + + if ( !fd.hasId(30) || !fd.hasId(31) || !fd.hasId(32) ) + pid_ploads = gsMatrix::Zero(1,points.cols()); + + for (index_t k =0; k!=points.cols(); k++) + pLoads.addLoad(points.col(k), loads.col(k), pid_ploads.at(k) ); // in parametric domain! + + + // Reference points + gsMatrix refPatches; + gsMatrix<> refPoints, refPars, refValue; // todo: add refValue.. + gsInfo<<"Reading reference point locations (ID=50) ..."; + if ( fd.hasId(50) ) fd.getId(50,refPoints); + gsInfo<<"Finished\n"; + gsInfo<<"Reading reference patches (ID=51) ..."; + if ( fd.hasId(51) ) fd.getId(51,refPatches); + gsInfo<<"Finished\n"; + gsInfo<<"Reading reference values (ID=52) ..."; + if ( fd.hasId(52) ) fd.getId(52,refValue); + gsInfo<<"Finished\n"; + + if ( !fd.hasId(50) || !fd.hasId(51) || !fd.hasId(52) ) + refValue = gsMatrix<>::Zero(mp.geoDim(),refPoints.cols()); + GISMO_ENSURE(refPatches.cols()==refPoints.cols(),"Number of reference points and patches do not match"); + + if (refPoints.rows()==2) + { + refPars = refPoints; + gsInfo<<"Reference points are provided in parametric coordinates.\n"; + } + else if (refPoints.rows()==3) + gsInfo<<"Reference points are provided in physical coordinates.\n"; + else + gsInfo<<"No reference points are provided.\n"; + + gsInfo<<"Reading DR solver options (ID=90) ..."; + if ( fd.hasId(90) ) fd.getId(90,DROptions); + gsInfo<<"Finished\n"; + gsInfo<<"Reading NR solver options (ID=91) ..."; + if ( fd.hasId(91) ) fd.getId(91,NROptions); + gsInfo<<"Finished\n"; + gsInfo<<"Reading assembler options (ID=92) ..."; + if ( fd.hasId(92) ) fd.getId(92,assemblerOptions); + gsInfo<<"Finished\n"; + + // ![Read data] + + // Fix path + dirname = gsFileManager::getCanonicRepresentation(dirname,true); + char sep = gsFileManager::getNativePathSeparator(); + gsFileManager::mkdir(dirname); + + // plot geometry + if (plot) + gsWriteParaview(mp,"mp",1000,true); + + mp.computeTopology(); + + // Make unstructured spline + gsMultiPatch<> geom; + gsMappedBasis<2,real_t> bb2; + gsSparseMatrix<> global2local; + if (method==0) + { + gsSmoothInterfaces<2,real_t> smoothInterfaces(mp); + smoothInterfaces.options().setSwitch("SharpCorners",false); + smoothInterfaces.compute(); + smoothInterfaces.matrix_into(global2local); + + global2local = global2local.transpose(); + geom = smoothInterfaces.exportToPatches(); + dbasis = smoothInterfaces.localBasis(); + } + else if (method==1) + { + gsAlmostC1<2,real_t> almostC1(mp); + almostC1.options().setSwitch("SharpCorners",true); + almostC1.compute(); + almostC1.matrix_into(global2local); + + global2local = global2local.transpose(); + geom = almostC1.exportToPatches(); + dbasis = almostC1.localBasis(); + } + else if (method==2) + { + gsDPatch<2,real_t> dpatch(mp); + dpatch.options().setSwitch("SharpCorners",true); + dpatch.compute(); + dpatch.matrix_into(global2local); + + global2local = global2local.transpose(); + geom = dpatch.exportToPatches(); + dbasis = dpatch.localBasis(); + } + else + GISMO_ERROR("Method "<* assembler; + if (membrane) + assembler = new gsThinShellAssembler<3, real_t, false >(geom,dbasis,BCs,forceFun,materialMatrixContainer); + else + assembler = new gsThinShellAssembler<3, real_t, true >(geom,dbasis,BCs,forceFun,materialMatrixContainer); + assembler->setOptions(assemblerOptions); + assembler->options().setInt("Continuity",-1); + assembler->setSpaceBasis(bb2); + assembler->setPointLoads(pLoads); + if (pressure) + assembler->setPressure(pressFun); + + // Assemble linear system to obtain the force vector + gsDebugVar(assembler->numDofs()); + assembler->assemble(); + gsSparseMatrix<> K = assembler->matrix(); + gsVector<> F = assembler->rhs(); + assembler->assembleMass(true); + gsVector<> M = assembler->rhs(); + + gsStructuralAnalysisOps::Jacobian_t Jacobian; + gsStructuralAnalysisOps::Residual_t Residual; + // Function for the Jacobian + Jacobian = [&assembler,&bb2,&geom](gsVector const &x, gsSparseMatrix & m) + { + ThinShellAssemblerStatus status; + gsMatrix solFull = assembler->fullSolutionVector(x); + size_t d = geom.targetDim(); + GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); + solFull.resize(solFull.rows()/d,d); + gsMappedSpline<2,real_t> mspline(bb2,solFull); + gsFunctionSum def(&geom,&mspline); + assembler->assembleMatrix(def); + status = assembler->assembleMatrix(def); + m = assembler->matrix(); + return status == ThinShellAssemblerStatus::Success; + }; + // Function for the Residual + Residual = [&assembler,&bb2,&geom](gsVector const &x, gsVector & result) + { + ThinShellAssemblerStatus status; + gsMatrix solFull = assembler->fullSolutionVector(x); + size_t d = geom.targetDim(); + GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); + solFull.resize(solFull.rows()/d,d); + + gsMappedSpline<2,real_t> mspline(bb2,solFull); + gsFunctionSum def(&geom,&mspline); + + status = assembler->assembleVector(def); + result = assembler->rhs(); // assembler rhs - force = Finternal + return status == ThinShellAssemblerStatus::Success; + }; + + gsStructuralAnalysisOutput writer(dirname + sep + wn,refPoints); + std::vector pointheaders = {"x","y","z"}; + std::vector otherheaders = {}; + writer.init(pointheaders,otherheaders); + + gsStaticNewton LIN(K,F); + LIN.setOptions(NROptions); + + gsStaticDR DRM(M,F,Residual); + DRM.setOptions(DROptions); + + gsStaticNewton NRM(K,F,Jacobian,Residual); + NRM.setOptions(NROptions); + + std::vector *> solvers{&LIN}; + if (DR) + { + gsInfo<<"Using Dynamic Relaxation solver\n"; + solvers.push_back(&DRM); + } + if (NR) + { + gsInfo<<"Using Newton-Raphson solver\n"; + solvers.push_back(&NRM); + } + gsStaticComposite solver(solvers); + solver.initialize(); + solver.solve(); + GISMO_ASSERT(solver.converged(),"Solver failed"); + gsVector<> solVector = solver.solution(); + + // ! [Export visualization in ParaView] + if (plot) + { + gsInfo<<"Plotting in Paraview...\n"; + std::string fileName; + /// Make a gsMappedSpline to represent the solution + // 1. Get all the coefficients (including the ones from the eliminated BCs.) + gsMatrix solFull = assembler->fullSolutionVector(solVector); + + // 2. Reshape all the coefficients to a Nx3 matrix + size_t d = geom.targetDim(); + GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); + solFull.resize(solFull.rows()/d,d); + + // 3. Make the mapped spline + gsMappedSpline<2,real_t> mspline(bb2,solFull); + + // 4. Create deformation spline + gsFunctionSum def(&geom,&mspline); + + // 5. Plot the mapped spline on the original geometry + gsPiecewiseFunction<> displacements; + assembler->constructStress(def,displacements,stress_type::displacement); + fileName = dirname + sep + "solution"; + gsWriteParaview(def,displacements,fileName,1000,"_"); + + // 5. Construct stress + gsPiecewiseFunction<> tensionFields; + assembler->constructStress(def,tensionFields,stress_type::tension_field); + fileName = dirname + sep + "tensionfield"; + gsWriteParaview(def,tensionFields,fileName,1000,"_"); + + + // gsExprAssembler<> A(1,1); + // A.setIntegrationElements(dbasis); + // gsExprEvaluator<> evA(A); + // gsExprAssembler<>::space u = A.getSpace(bb2,3); + + // // gsConstantFunction<> one(1.0,3); + // // gsBoundaryConditions<> bc_dummy; + // // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + + // gsFunctionExpr<> Mz_x("sqrt(x^2+y^2)*cos(atan2(y,x))",3); + // gsFunctionExpr<> Mz_y("sqrt(x^2+y^2)*sin(atan2(y,x))",3); + // gsBoundaryConditions<> bc_dummy; + // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + + // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + + + // bc_dummy.setGeoMap(geom); + // u.setup(bc_dummy,dirichlet::interpolation,-1); + // gsDebugVar(u.mapper()); + + // gsMatrix<> coefs(bb2.size(),3); + // for (size_t d=0; d!=geom.geoDim(); d++) + // for (index_t k=0; k!=bb2.size(); k++) + // { + // const index_t ii = u.mapper().index(k,0,d); + // if (u.mapper().is_free(k,0,d)) + // coefs(k,d) = 0.0; + // else if (u.mapper().is_boundary(k,0,d)) + // coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); + // else + // GISMO_ERROR("WHAT?"); + // } + + // gsMappedSpline<2,real_t> geo(bb2,coefs); + + + gsExprAssembler<> A(1,1); + A.setIntegrationElements(dbasis); + gsExprEvaluator<> evA(A); + gsExprAssembler<>::space u = A.getSpace(dbasis,3); + + // gsConstantFunction<> one(1.0,3); // for X reaction force + // gsBoundaryConditions<> bc_dummy; + // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + + gsBoundaryConditions<> bc_dummy; + gsFunctionExpr<> Mz_x("sqrt(x^2+y^2)*cos(atan2(y,x))",3); + gsFunctionExpr<> Mz_y("-sqrt(x^2+y^2)*sin(atan2(y,x))",3); // X loads contribute negatively to M_z, so we need to flip the sign + bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + + bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + +/* + gsBoundaryConditions<> bc_dummy; + gsMultiPatch<> mp_x = mp.coord(0); + gsMultiPatch<> mp_y = mp.coord(1); + // X loads contribute negatively to M_z, so we need to flip the sign + for (size_t p=0; p!=mp_y.nPatches(); p++) + mp_y.patch(p).coefs().array() *= -1; + bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); + bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); + bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); + bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); + + bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); + bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); + bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); + bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); +*/ + bc_dummy.setGeoMap(mp); + u.setup(bc_dummy,dirichlet::interpolation,-1); + gsDebugVar(u.mapper()); + + gsMultiPatch<> geo; + for (size_t b=0; b!=dbasis.nBases(); b++) + { + gsMatrix<> coefs(dbasis.basis(b).size(),mp.geoDim()); + coefs.setZero(); + for (size_t d=0; d!=mp.geoDim(); d++) + for (index_t k=0; k!=dbasis.basis(b).size(); k++) + { + const index_t ii = u.mapper().index(k,b,d); + if (u.mapper().is_free(k,b,d)) + coefs(k,d) = 0.0; + else if (u.mapper().is_boundary(k,b,d)) + coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); + else + GISMO_ERROR("WHAT?"); + } + geo.addPatch(*dbasis.basis(b).makeGeometry(give(coefs))); + + } + + // for (size_t b=0; b!=dbasis.nBases(); b++) + // for (size_t d=0; d!=2; d++) + // for (index_t k=0; k!=dbasis.basis(b).size(); k++) + // { + // const index_t ii = u.mapper().index(k,0,d); + // if (u.mapper().is_free(k,0,d)) + // coefs(k,d) = 0.0; + // else if (u.mapper().is_boundary(k,0,d)) + // coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); + // else + // GISMO_ERROR("WHAT?"); + // } + + // gsGeometry<>::uPtr geo = dbasis.basis(0).makeGeometry(give(coefs)); + + A.initSystem(); + + gsExprAssembler<>::geometryMap ori = A.getMap(mp); + gsExprAssembler<>::geometryMap deff = A.getMap(def); + + auto v = A.getCoeff(geo); + + gsMaterialMatrixIntegrate S0f(materialMatrixContainer,&mp,&def); + gsMaterialMatrixIntegrate S1f(materialMatrixContainer,&mp,&def); + auto S0 = A.getCoeff(S0f); + auto S1 = A.getCoeff(S1f); + + // Helper matrix for flexural components + gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2); + auto m2 = A.getCoeff(mult2t); + + auto N = S0.tr(); + auto dEm= flat( jac(deff).tr() * jac(v) ) ; + auto M = S1.tr(); + auto dEf= -( deriv2(v,sn(deff).normalized().tr() ) + deriv2(deff,var1(v,deff) ) ) * reshape(m2,3,3); + + + gsDebugVar(evA.integral(( + N*dEm.tr() + + M*dEf.tr() + )*meas(ori))); + // std::vector bContainer({patchSide(0,boundary::north),patchSide(1,boundary::north),patchSide(2,boundary::north),patchSide(3,boundary::north)}); + + + + + + + } + if (write) + { + // 1. Get all the coefficients (including the ones from the eliminated BCs.) + gsMatrix solFull = assembler->fullSolutionVector(solVector); + + // 2. Reshape all the coefficients to a Nx3 matrix + size_t d = geom.targetDim(); + GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); + solFull.resize(solFull.rows()/d,d); + + // 3. Make the mapped spline + gsMappedSpline<2,real_t> mspline(bb2,solFull); + + // 4. Make geom_def by projecting coefficients + gsDofMapper mapper(dbasis); + mapper.finalize(); + gsMatrix<> coefs; + gsMultiPatch<> geom_def = geom; + gsInfo<<"L2-Projection error of mspline on dbasis = "<::projectFunction(dbasis,mspline,geom,coefs)<<"\n"; + coefs.resize(coefs.rows()/geom.geoDim(),geom.geoDim()); + + index_t offset = 0; + for (size_t p = 0; p != geom_def.nPatches(); p++) + { + gsMatrix<> tmp_coefs = geom.patch(p).coefs(); + tmp_coefs += coefs.block(offset,0,mapper.patchSize(p),geom.geoDim()); + geom_def.patch(p) = give(*dbasis.basis(p).makeGeometry((tmp_coefs))); + offset += mapper.patchSize(p); + } + + // 5. Write the result + gsWrite(geom_def,dirname + sep + "deformed"); + + // gsMatrix<> result(mp.geoDim(),refPoints.cols()); + // for (index_t k=0; k!=refPoints.cols(); k++) + // result.col(k) = deformation.patch(refPatches.at(k)).eval(refPoints.col(k)); + // writer.add(result); + } + + // ! [Export visualization in ParaView] + if (stress) + { + gsField<> tensionField; + gsPiecewiseFunction<> tensionFields; + std::string fileName = dirname + sep + "tensionfield"; + + /// Make a gsMappedSpline to represent the solution + // 1. Get all the coefficients (including the ones from the eliminated BCs.) + gsMatrix solFull = assembler->fullSolutionVector(solVector); + + // 2. Reshape all the coefficients to a Nx3 matrix + size_t d = geom.targetDim(); + GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); + solFull.resize(solFull.rows()/d,d); + + // 3. Make the mapped spline + gsMappedSpline<2,real_t> mspline(bb2,solFull); + + // 4. Create deformation spline + gsFunctionSum def(&geom,&mspline); + + // 5. Construct stress + assembler->constructStress(def,tensionFields,stress_type::tension_field); + gsWriteParaview(def,tensionFields,fileName,1000,"_"); + + + // gsPiecewiseFunction<> membraneStresses; + // assembler->constructStress(mp_def,membraneStresses,stress_type::membrane); + // gsField<> membraneStress(mp_def,membraneStresses, true); + + // gsPiecewiseFunction<> flexuralStresses; + // assembler->constructStress(mp_def,flexuralStresses,stress_type::flexural); + // gsField<> flexuralStress(mp_def,flexuralStresses, true); + + // gsPiecewiseFunction<> stretches; + // assembler->constructStress(mp_def,stretches,stress_type::principal_stretch); + // gsField<> Stretches(mp_def,stretches, true); + + // gsPiecewiseFunction<> pstrain_m; + // assembler->constructStress(mp_def,pstrain_m,stress_type::principal_membrane_strain); + // gsField<> pstrainM(mp_def,pstrain_m, true); + + // gsPiecewiseFunction<> pstrain_f; + // assembler->constructStress(mp_def,pstrain_f,stress_type::principal_flexural_strain); + // gsField<> pstrainF(mp_def,pstrain_f, true); + + // gsPiecewiseFunction<> pstress_m; + // assembler->constructStress(mp_def,pstress_m,stress_type::principal_stress_membrane); + // gsField<> pstressM(mp_def,pstress_m, true); + + // gsPiecewiseFunction<> pstress_f; + // assembler->constructStress(mp_def,pstress_f,stress_type::principal_stress_flexural); + // gsField<> pstressF(mp_def,pstress_f, true); + + // gsPiecewiseFunction<> stretch1; + // assembler->constructStress(mp_def,stretch1,stress_type::principal_stretch_dir1); + // gsField<> stretchDir1(mp_def,stretch1, true); + + // gsPiecewiseFunction<> stretch2; + // assembler->constructStress(mp_def,stretch2,stress_type::principal_stretch_dir2); + // gsField<> stretchDir2(mp_def,stretch2, true); + + // gsPiecewiseFunction<> stretch3; + // assembler->constructStress(mp_def,stretch3,stress_type::principal_stretch_dir3); + // gsField<> stretchDir3(mp_def,stretch3, true); + + // gsPiecewiseFunction<> VMStresses; + // assembler->constructStress(mp_def,VMStresses,stress_type::von_mises_membrane); + // gsField<> VMStress(mp_def,VMStresses, true); + + + // gsWriteParaview(membraneStress,dirname + sep + "MembraneStress",5000); + // gsWriteParaview(VMStress,dirname + sep + "MembraneStressVM",5000); + // gsWriteParaview(Stretches,dirname + sep + "PrincipalStretch",5000); + // gsWriteParaview(pstrainM,dirname + sep + "PrincipalMembraneStrain",5000); + // gsWriteParaview(pstrainF,dirname + sep + "PrincipalFlexuralStrain",5000); + // gsWriteParaview(pstressM,dirname + sep + "PrincipalMembraneStress",5000); + // gsWriteParaview(pstressF,dirname + sep + "PrincipalFlexuralStress",5000); + // gsWriteParaview(stretchDir1,dirname + sep + "PrincipalDirection1",5000); + // gsWriteParaview(stretchDir1,dirname + sep + "PrincipalDirection1",5000); + // gsWriteParaview(stretchDir2,dirname + sep + "PrincipalDirection2",5000); + // gsWriteParaview(stretchDir3,dirname + sep + "PrincipalDirection3",5000); + + } + + delete assembler; + + return EXIT_SUCCESS; +} +#else//gsUnstructuredSplines_ENABLED +int main(int argc, char *argv[]) +{ + gsWarn<<"G+Smo is not compiled with the gsUnstructuredSplines module."; + return EXIT_FAILURE; +} +#endif +#else//gsKLShell_ENABLED +int main(int argc, char *argv[]) +{ + gsWarn<<"G+Smo is not compiled with the gsKLShell module."; + return EXIT_FAILURE; +} +#endif \ No newline at end of file From 2f144667b55d30b6ab13ce777be2b2dc33b9e084 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Fri, 29 Nov 2024 18:42:06 +0100 Subject: [PATCH 02/29] add explicit case for cylinder and annulus --- benchmarks/benchmark_cylinder_DC.cpp | 552 +++++++++++++++++++++++++++ src/gsStaticSolvers/gsStaticNewton.h | 6 +- 2 files changed, 557 insertions(+), 1 deletion(-) create mode 100644 benchmarks/benchmark_cylinder_DC.cpp diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp new file mode 100644 index 0000000..4a61c82 --- /dev/null +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -0,0 +1,552 @@ + +/** @file benchmark_TensionWrinkling.cpp + + @brief Computes the wrinkling behaviour of a thin sheet + + 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 (2019-..., TU Delft) +*/ + +#include + +#ifdef gsKLShell_ENABLED +#include +#include +#include +#include +#endif + +#ifdef gsUnstructuredSplines_ENABLED +#include +#endif + +#include +#include +#include +#include +using namespace gismo; + +#ifdef gsKLShell_ENABLED + +int main (int argc, char** argv) +{ + // Input options + int numElevate = 2; + int numHref = 5; + bool plot = false; + bool mesh = false; + int step = 10; + + bool TFT = false; + real_t dL = 1; + // Solvers + // * Dynamic Relaxation + index_t maxitDR = 1e4; + real_t alpha = 1.0; + real_t damping = 0.1; + bool DR = false; + // * Newton Raphson + index_t maxitNR = 20; + bool NR = false; + index_t testCase = 0; + + + std::string wn("data.csv"); + + std::string assemberOptionsFile("options/solver_options.xml"); + + gsCmdLine cmd("Wrinkling analysis with thin shells."); + // Input settings + cmd.addString( "f", "file", "Input XML file for assembler options", assemberOptionsFile ); + + // Mesh settings + cmd.addInt("r","hRefine", "Number of dyadic h-refinement (bisection) steps to perform before solving", numHref); + cmd.addInt("e","degreeElevation", "Number of degree elevation steps to perform on the Geometry's basis before solving", numElevate); + + // Material settings + cmd.addSwitch("TFT", "Use TFT material matrix", TFT); + + // Test case + cmd.addInt("t", "test", "Test case: (0): Annulus | (1): Cylinder", testCase); + + // Load/displ stepping settings + cmd.addInt("N", "maxsteps", "Maximum number of steps", step); + cmd.addReal( "L", "step", "load step size", dL); + cmd.addReal( "a", "alpha", "alpha", alpha ); + cmd.addReal( "c", "damping", "damping", damping ); + cmd.addSwitch("DR", "Use Dynamic Relaxation", DR); + cmd.addSwitch("NR", "Use Newton Raphson", NR); + + // Output settings + cmd.addSwitch("plot", "Plot result in ParaView format", plot); + cmd.addSwitch("mesh", "Plot mesh?", mesh); + + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + std::string output = "solution"; + std::string dirname = "ArcLengthResults"; + + gsFileData<> fd(assemberOptionsFile); + gsOptionList opts; + fd.getFirst(opts); + + real_t R = 0.25, H = 1.0; + real_t thickness = 0.05e-3; + real_t YoungsModulus = 1e9; + real_t PoissonRatio = 0.5; + real_t Density = 1e0; + + gsMultiPatch<> mp, mp_def; + if (testCase==0) + { + mp = gsNurbsCreator<>::MultiPatchAnnulus(0.0625,0.25); + mp.embed(3); + } + else if (testCase==1) + mp = gsNurbsCreator<>::MultiPatchCylinder(0.25,1.0); + + mp.computeTopology(); + + for (size_t p = 0; p!=mp.nPatches(); ++p) + { + for(index_t i = 0; i< numElevate; ++i) + mp.patch(p).degreeElevate(); // Elevate the degree + for(index_t i = 0; i< numHref; ++i) + mp.patch(p).uniformRefine(); + } + + mp_def = mp; + + // Make unstructured spline + gsMultiPatch<> geom; + gsMultiBasis<> dbasis; + gsMappedBasis<2,real_t> bb2; + gsSparseMatrix<> global2local; + gsSmoothInterfaces<2,real_t> smoothInterfaces(mp); + smoothInterfaces.options().setSwitch("SharpCorners",false); + smoothInterfaces.compute(); + smoothInterfaces.matrix_into(global2local); + + global2local = global2local.transpose(); + geom = smoothInterfaces.exportToPatches(); + dbasis = smoothInterfaces.localBasis(); + bb2.init(dbasis,global2local); + + // Assign BCs + gsFunctionExpr<> dx("sqrt(x^2+y^2)*(cos(atan2(y,x) + pi/2*t)- cos(atan2(y,x)))",3); + gsFunctionExpr<> dy("sqrt(x^2+y^2)*(sin(atan2(y,x) + pi/2*t)- sin(atan2(y,x)))",3); + real_t dz_val = 0; + if (testCase==0) + dz_val = 125e-3; + else if (testCase==1) + dz_val = 1.0; + + gsConstantFunction<> dz(dz_val,3); + + // Boundary conditions + gsBoundaryConditions<> BCs; + BCs.setGeoMap(geom); + for (size_t p = 0; p!=geom.nPatches(); ++p) + { + BCs.addCondition(p,boundary::south, condition_type::dirichlet, 0 , 0, false, -1); + BCs.addCondition(p,boundary::north, condition_type::dirichlet, &dx, 0, false, 0); + BCs.addCondition(p,boundary::north, condition_type::dirichlet, &dy, 0, false, 1); + BCs.addCondition(p,boundary::north, condition_type::dirichlet, &dz, 0, false, 2); + } + + if (testCase==0) + dirname = dirname + "/Annulus_-r" + std::to_string(numHref) + "-e" + std::to_string(numElevate); + else if (testCase==1) + dirname = dirname + "/Cylinder_-r" + std::to_string(numHref) + "-e" + std::to_string(numElevate); + else + GISMO_ERROR("Test case "< force("0","0","0",3); + + gsConstantFunction<> t(thickness,3); + gsConstantFunction<> E(YoungsModulus,3); + gsConstantFunction<> nu(PoissonRatio,3); + gsConstantFunction<> rho(Density,3); + + std::vector*> parameters; + parameters.resize(2); + parameters[0] = &E; + parameters[1] = ν + + gsMaterialMatrixBase * materialMatrix; + gsMaterialMatrixBase * materialMatrixTFT; + gsThinShellAssemblerBase* assembler; + + gsOptionList options; + options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",1); + options.addSwitch("Compressibility","Compressibility: (false): Imcompressible | (true): Compressible",false); + options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",1); + materialMatrix = getMaterialMatrix<3,real_t>(mp,t,parameters,rho,options); + gsMaterialMatrixContainer materialMatrixContainer; + if (TFT) + { + materialMatrixTFT = new gsMaterialMatrixTFT<3,real_t,true>(materialMatrix); + for (size_t p = 0; p!=mp.nPatches(); p++) + materialMatrixContainer.add(materialMatrixTFT); + // materialMatrixTFT->options().setReal("SlackMultiplier",1e-6); + assembler = new gsThinShellAssembler<3, real_t, false >(geom,dbasis,BCs,force,materialMatrixContainer); + } + else + { + for (size_t p = 0; p!=mp.nPatches(); p++) + materialMatrixContainer.add(materialMatrix); + assembler = new gsThinShellAssembler<3, real_t, true >(geom,dbasis,BCs,force,materialMatrixContainer); + } + + + assembler->setOptions(opts); + assembler->options().setInt("Continuity",-1); + assembler->setSpaceBasis(bb2); + + // Assemble linear system to obtain the force vector + assembler->assemble(); + gsVector<> F = assembler->rhs(); + gsSparseMatrix<> K = assembler->matrix(); + assembler->assembleMass(true); + gsVector<> M = assembler->rhs(); + + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&bb2, &geom](gsVector const &x, gsSparseMatrix & m) + { + ThinShellAssemblerStatus status; + gsMatrix solFull = assembler->fullSolutionVector(x); + size_t d = geom.targetDim(); + GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); + solFull.resize(solFull.rows()/d,d); + gsMappedSpline<2,real_t> mspline(bb2,solFull); + gsFunctionSum def(&geom,&mspline); + assembler->assembleMatrix(def); + status = assembler->assembleMatrix(def); + m = assembler->matrix(); + return status == ThinShellAssemblerStatus::Success; + }; + + // Function for the Residual + gsStructuralAnalysisOps::Residual_t Residual = [&dx,&dy,&BCs,&assembler,&bb2, &geom](gsVector const &x, gsVector & result) + { + ThinShellAssemblerStatus status; + gsMatrix solFull = assembler->fullSolutionVector(x); + size_t d = geom.targetDim(); + GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); + solFull.resize(solFull.rows()/d,d); + + gsMappedSpline<2,real_t> mspline(bb2,solFull); + gsFunctionSum def(&geom,&mspline); + + status = assembler->assembleVector(def); + result = assembler->rhs(); + return status == ThinShellAssemblerStatus::Success; + + }; + + gsParaviewCollection collection(dirname + "/" + output); + gsParaviewCollection collectionTF(dirname + "/" + "tensionfield"); + gsMultiPatch<> deformation = mp; + + // Define solvers + gsStaticDR DynamicRelaxationSolver(M,F,Residual); + gsOptionList DynamicRelaxationOptions = DynamicRelaxationSolver.options(); + DynamicRelaxationOptions.setReal("damping",damping); + DynamicRelaxationOptions.setReal("alpha",alpha); + DynamicRelaxationOptions.setInt("maxIt",maxitDR); + DynamicRelaxationOptions.setReal("tolE",1e-4); + DynamicRelaxationOptions.setReal("tol",1e-2); + DynamicRelaxationOptions.setInt("verbose",1); + DynamicRelaxationSolver.setOptions(DynamicRelaxationOptions); + DynamicRelaxationSolver.initialize(); + + gsStaticNewton NewtonSolver(K,F,Jacobian,Residual); + gsOptionList NewtonOptions = NewtonSolver.options(); + NewtonOptions.setInt("verbose",true); + NewtonOptions.setInt("maxIt",maxitNR); + NewtonOptions.setReal("tolU",1e-4); + NewtonOptions.setReal("tolF",1e-6); + NewtonSolver.setOptions(NewtonOptions); + + std::vector *> solvers; + if (DR) + { + gsInfo<<"Using Dynamic Relaxation solver\n"; + solvers.push_back(&DynamicRelaxationSolver); + } + if (NR) + { + gsInfo<<"Using Newton-Raphson solver\n"; + solvers.push_back(&NewtonSolver); + } + + gsStaticComposite StaticSolver(solvers); + // gsControlDisplacement controlDC(&StaticSolver); + // Displacement-controlled simulation + real_t dL0 = dL; + real_t D = dL; + gsMatrix<> U = gsMatrix<>::Zero(assembler->numDofs(),1); + + std::string writeName = dirname + "/out.csv"; + std::ofstream file; + file.open(writeName); + file<<"Load factor, Moment\n"; + file<updateBCs(BCs); + StaticSolver.setDisplacement(U); + StaticSolver.solve(); + U = StaticSolver.solution(); + + // 1. Get all the coefficients (including the ones from the eliminated BCs.) + // gsMatrix solFull = assembler->fullSolutionVector(controlDC.solutionU()); + gsMatrix solFull = assembler->fullSolutionVector(U); + + // 2. Reshape all the coefficients to a Nx3 matrix + size_t d = geom.targetDim(); + GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); + solFull.resize(solFull.rows()/d,d); + + // 3. Make the mapped spline + gsMappedSpline<2,real_t> mspline(bb2,solFull); + + // 4. Create deformation spline + gsFunctionSum def(&geom,&mspline); + + if (plot) + { + // 5. Plot the mapped spline on the original geometry + gsPiecewiseFunction<> displacements; + assembler->constructStress(def,displacements,stress_type::displacement); + std::string fileName = dirname + "/" + output + util::to_string(k); + gsWriteParaview(def,displacements,fileName,10000,"_"); + fileName = gsFileManager::getFilename(fileName); + for (size_t p=0; p tensionFields; + assembler->constructStress(def,tensionFields,stress_type::tension_field); + std::string fileNameTF = dirname + "/" + "tensionfield" + util::to_string(k); + gsWriteParaview(def,tensionFields,fileNameTF,10000,"_"); + fileNameTF = gsFileManager::getFilename(fileNameTF); + for (size_t p=0; p A(1,1); + // A.setIntegrationElements(dbasis); + // gsExprEvaluator<> evA(A); + // gsExprAssembler<>::space u = A.getSpace(bb2,3); + + // // gsConstantFunction<> one(1.0,3); + // // gsBoundaryConditions<> bc_dummy; + // // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + + // gsFunctionExpr<> Mz_x("sqrt(x^2+y^2)*cos(atan2(y,x))",3); + // gsFunctionExpr<> Mz_y("sqrt(x^2+y^2)*sin(atan2(y,x))",3); + // gsBoundaryConditions<> bc_dummy; + // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + + // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + + + // bc_dummy.setGeoMap(geom); + // u.setup(bc_dummy,dirichlet::interpolation,-1); + // gsDebugVar(u.mapper()); + + // gsMatrix<> coefs(bb2.size(),3); + // for (size_t d=0; d!=geom.geoDim(); d++) + // for (index_t k=0; k!=bb2.size(); k++) + // { + // const index_t ii = u.mapper().index(k,0,d); + // if (u.mapper().is_free(k,0,d)) + // coefs(k,d) = 0.0; + // else if (u.mapper().is_boundary(k,0,d)) + // coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); + // else + // GISMO_ERROR("WHAT?"); + // } + + // gsMappedSpline<2,real_t> geo(bb2,coefs); + + + gsExprAssembler<> A(1,1); + A.setIntegrationElements(dbasis); + gsExprEvaluator<> evA(A); + gsExprAssembler<>::space u = A.getSpace(dbasis,3); + + // gsConstantFunction<> one(1.0,3); // for X reaction force + // gsBoundaryConditions<> bc_dummy; + // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); + + gsBoundaryConditions<> bc_dummy; + gsFunctionExpr<> Mz_x("sqrt(x^2+y^2)*cos(atan2(y,x))",3); + gsFunctionExpr<> Mz_y("-sqrt(x^2+y^2)*sin(atan2(y,x))",3); // X loads contribute negatively to M_z, so we need to flip the sign + bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); + + bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); + + /* + gsBoundaryConditions<> bc_dummy; + gsMultiPatch<> mp_x = mp.coord(0); + gsMultiPatch<> mp_y = mp.coord(1); + // X loads contribute negatively to M_z, so we need to flip the sign + for (size_t p=0; p!=mp_y.nPatches(); p++) + mp_y.patch(p).coefs().array() *= -1; + bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); + bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); + bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); + bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); + + bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); + bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); + bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); + bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); + */ + bc_dummy.setGeoMap(mp); + u.setup(bc_dummy,dirichlet::interpolation,-1); + gsDebugVar(u.mapper()); + + gsMultiPatch<> geo; + for (size_t b=0; b!=dbasis.nBases(); b++) + { + gsMatrix<> coefs(dbasis.basis(b).size(),mp.geoDim()); + coefs.setZero(); + for (size_t d=0; d!=mp.geoDim(); d++) + for (index_t k=0; k!=dbasis.basis(b).size(); k++) + { + const index_t ii = u.mapper().index(k,b,d); + if (u.mapper().is_free(k,b,d)) + coefs(k,d) = 0.0; + else if (u.mapper().is_boundary(k,b,d)) + coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); + else + GISMO_ERROR("WHAT?"); + } + geo.addPatch(*dbasis.basis(b).makeGeometry(give(coefs))); + + } + + // for (size_t b=0; b!=dbasis.nBases(); b++) + // for (size_t d=0; d!=2; d++) + // for (index_t k=0; k!=dbasis.basis(b).size(); k++) + // { + // const index_t ii = u.mapper().index(k,0,d); + // if (u.mapper().is_free(k,0,d)) + // coefs(k,d) = 0.0; + // else if (u.mapper().is_boundary(k,0,d)) + // coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); + // else + // GISMO_ERROR("WHAT?"); + // } + + // gsGeometry<>::uPtr geo = dbasis.basis(0).makeGeometry(give(coefs)); + + A.initSystem(); + + gsExprAssembler<>::geometryMap ori = A.getMap(mp); + gsExprAssembler<>::geometryMap deff = A.getMap(def); + + auto v = A.getCoeff(geo); + + gsMaterialMatrixIntegrate S0f(materialMatrixContainer,&mp,&def); + gsMaterialMatrixIntegrate S1f(materialMatrixContainer,&mp,&def); + auto S0 = A.getCoeff(S0f); + auto S1 = A.getCoeff(S1f); + + // Helper matrix for flexural components + gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2); + auto m2 = A.getCoeff(mult2t); + + auto N = S0.tr(); + auto dEm= flat( jac(deff).tr() * jac(v) ) ; + auto M = S1.tr(); + auto dEf= -( deriv2(v,sn(deff).normalized().tr() ) + deriv2(deff,var1(v,deff) ) ) * reshape(m2,3,3); + + real_t Moment = evA.integral(( + N*dEm.tr() + + M*dEf.tr() + )*meas(ori) + ); + + file.open(writeName,std::ofstream::out | std::ofstream::app); + file<< std::setprecision(10) + << D << "," + << Moment << "\n"; + file.close(); + + + + + // controlDC.step(dL); + D+= dL; + } + + if (plot) + { + collection.save(); + collectionTF.save(); + } + + delete materialMatrix; + delete assembler; + if (TFT) + delete materialMatrixTFT; + return EXIT_SUCCESS; +} + + + +#else//gsKLShell_ENABLED +int main(int argc, char *argv[]) +{ + gsWarn<<"G+Smo is not compiled with the gsKLShell module."; + return EXIT_FAILURE; +} +#endif \ No newline at end of file diff --git a/src/gsStaticSolvers/gsStaticNewton.h b/src/gsStaticSolvers/gsStaticNewton.h index 80ad71f..6263796 100644 --- a/src/gsStaticSolvers/gsStaticNewton.h +++ b/src/gsStaticSolvers/gsStaticNewton.h @@ -55,6 +55,7 @@ class gsStaticNewton : public gsStaticBase m_nonlinear(nullptr), m_residualFun(nullptr) { + m_dofs = m_linear.rows(); this->_init(); } @@ -78,6 +79,7 @@ class gsStaticNewton : public gsStaticBase m_residualFun(residual), m_ALresidualFun(nullptr) { + m_dofs = m_linear.rows(); m_dnonlinear = [this](gsVector const & x, gsVector const & dx, gsSparseMatrix & m) -> bool { return m_nonlinear(x,m); @@ -106,6 +108,7 @@ class gsStaticNewton : public gsStaticBase m_residualFun(nullptr), m_ALresidualFun(ALResidual) { + m_dofs = m_linear.rows(); m_L = 1.0; m_residualFun = [this](gsVector const & x, gsVector & result) -> bool { @@ -141,6 +144,7 @@ class gsStaticNewton : public gsStaticBase m_residualFun(residual), m_ALresidualFun(nullptr) { + m_dofs = m_linear.rows(); this->_init(); } @@ -270,7 +274,7 @@ class gsStaticNewton : public gsStaticBase /// Linear solver employed using Base::m_solver; // Cholesky by default - + using Base::m_stabilityMethod; /// Indicator for bifurcation From 4c56065b3d6a533826f5d18f6eeec753496bdb12 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 4 Dec 2024 11:10:45 +0100 Subject: [PATCH 03/29] static solver changes --- src/gsStaticSolvers/gsStaticBase_.cpp | 16 ++++++++++++++++ src/gsStaticSolvers/gsStaticOpt.h | 8 ++++---- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/src/gsStaticSolvers/gsStaticBase_.cpp b/src/gsStaticSolvers/gsStaticBase_.cpp index d153fe3..67759c8 100644 --- a/src/gsStaticSolvers/gsStaticBase_.cpp +++ b/src/gsStaticSolvers/gsStaticBase_.cpp @@ -16,6 +16,9 @@ #ifdef gsHLBFGS_ENABLED #include #endif +#ifdef gsOptim_ENABLED +#include +#endif #include #include @@ -32,5 +35,18 @@ namespace gismo #ifdef gsHLBFGS_ENABLED CLASS_TEMPLATE_INST gsStaticOpt>; #endif +#ifdef gsOptim_ENABLED + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + // CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; +#endif } diff --git a/src/gsStaticSolvers/gsStaticOpt.h b/src/gsStaticSolvers/gsStaticOpt.h index 7800905..dd224a9 100644 --- a/src/gsStaticSolvers/gsStaticOpt.h +++ b/src/gsStaticSolvers/gsStaticOpt.h @@ -123,8 +123,8 @@ class gsStaticOpt : public gsStaticBase */ gsStaticOpt( const Residual_t &Residual, const index_t numDofs ) : - m_optProblem(Residual,m_L,numDofs), - m_optimizer(&m_optProblem) + m_optimizer(&m_optProblem), + m_optProblem(Residual,m_L,numDofs) { m_dofs = numDofs; this->_init(); @@ -137,8 +137,8 @@ class gsStaticOpt : public gsStaticBase */ gsStaticOpt( const ALResidual_t & ALResidual, const index_t numDofs ) : - m_optProblem(ALResidual,m_L,numDofs), - m_optimizer(&m_optProblem) + m_optimizer(&m_optProblem), + m_optProblem(ALResidual,m_L,numDofs) { m_L = 1.0; m_dofs = numDofs; From 338242c85e0cb30dbf5ea78bc8a6c38b8601ffba Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Thu, 5 Dec 2024 16:19:39 +0100 Subject: [PATCH 04/29] add `gsStaticOpt` in example only add solution in `gsStaticComposite` if converged --- benchmarks/benchmark_cylinder_DC.cpp | 30 +++++++++++++++++++++++ src/gsStaticSolvers/gsStaticComposite.hpp | 13 +++++++--- 2 files changed, 39 insertions(+), 4 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 4a61c82..2fb4d10 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -25,7 +25,16 @@ #include #endif +#ifdef gsHLBFGS_ENABLED +#include +#endif + +#ifdef gsOptim_ENABLED +#include +#endif + #include +#include #include #include #include @@ -50,9 +59,14 @@ int main (int argc, char** argv) real_t alpha = 1.0; real_t damping = 0.1; bool DR = false; + // * Optimizer + index_t maxitOPT = 1000; + bool OP = false; // * Newton Raphson index_t maxitNR = 20; bool NR = false; + + index_t testCase = 0; @@ -80,6 +94,7 @@ int main (int argc, char** argv) cmd.addReal( "a", "alpha", "alpha", alpha ); cmd.addReal( "c", "damping", "damping", damping ); cmd.addSwitch("DR", "Use Dynamic Relaxation", DR); + cmd.addSwitch("OP", "Use Optimizer", OP); cmd.addSwitch("NR", "Use Newton Raphson", NR); // Output settings @@ -277,6 +292,16 @@ int main (int argc, char** argv) DynamicRelaxationSolver.setOptions(DynamicRelaxationOptions); DynamicRelaxationSolver.initialize(); + gsStaticOpt::LBFGS> OptimizerSolver(Residual,assembler->numDofs()); + gsOptionList OptimizerSolverOptions = OptimizerSolver.options(); + OptimizerSolverOptions.setInt("maxIt",maxitOPT); + OptimizerSolverOptions.setReal("tolF",1e-4); + OptimizerSolverOptions.setReal("tolU",1e-4); + OptimizerSolverOptions.setInt("verbose",1); + OptimizerSolver.setOptions(OptimizerSolverOptions); + OptimizerSolver.initialize(); + + gsStaticNewton NewtonSolver(K,F,Jacobian,Residual); gsOptionList NewtonOptions = NewtonSolver.options(); NewtonOptions.setInt("verbose",true); @@ -291,6 +316,11 @@ int main (int argc, char** argv) gsInfo<<"Using Dynamic Relaxation solver\n"; solvers.push_back(&DynamicRelaxationSolver); } + if (OP) + { + gsInfo<<"Using Optimizer solver\n"; + solvers.push_back(&OptimizerSolver); + } if (NR) { gsInfo<<"Using Newton-Raphson solver\n"; diff --git a/src/gsStaticSolvers/gsStaticComposite.hpp b/src/gsStaticSolvers/gsStaticComposite.hpp index b84c3f6..61b7785 100644 --- a/src/gsStaticSolvers/gsStaticComposite.hpp +++ b/src/gsStaticSolvers/gsStaticComposite.hpp @@ -27,16 +27,21 @@ gsStatus gsStaticComposite::solve() { if (m_verbose > 0) gsInfo<<"Solver "<0) { m_solvers[k]->setUpdate(m_solvers[k-1]->update()); } m_status = m_solvers[k]->solve(); - m_numIterations += m_solvers[k]->iterations(); - m_U = m_solvers[k]->solution(); - m_DeltaU = m_solvers[k]->update(); + if (m_status==gsStatus::Success) + { + m_numIterations += m_solvers[k]->iterations(); + m_U = m_solvers[k]->solution(); + m_DeltaU = m_solvers[k]->update(); + } + else if (m_verbose > 0) + gsInfo<<"Solver "< Date: Sun, 8 Dec 2024 19:36:38 +0100 Subject: [PATCH 05/29] add unittest add docs static optimizer change cmakelists targets --- .github/workflows/ci.yml | 2 +- CMakeLists.txt | 8 +- benchmarks/CMakeLists.txt | 2 +- examples/CMakeLists.txt | 2 +- solvers/CMakeLists.txt | 2 +- src/gsStaticSolvers/gsStaticOpt.h | 88 ++++-- src/gsStaticSolvers/gsStaticOpt.hpp | 141 +--------- tutorials/CMakeLists.txt | 2 +- unittests/gsStaticSolver_test.cpp | 421 ++++++++++++++++++++++++++++ 9 files changed, 494 insertions(+), 174 deletions(-) create mode 100644 unittests/gsStaticSolver_test.cpp diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 21e42de..772260b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -49,4 +49,4 @@ jobs: - name: "Run for ${{ matrix.os }}" shell: bash working-directory: ${{runner.workspace}} - run: ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" -D CTEST_SITE="${{ matrix.os }}_[actions]" -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON" -D GISMO_OPTIONAL='gsSpectra\\;gsKLShell\\;gsElasticity\\;${{ github.event.repository.name }}' -Q + run: ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" -D CTEST_CONFIGURATION_TYPE=RelWithDebInfo -D LABELS_FOR_SUBPROJECTS="gsStructuralAnalysis-benchmarks" -D CTEST_SITE="${{ matrix.os }}_[actions]" -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON" -D GISMO_OPTIONAL='gsOptim\\;gsSpectra\\;gsKLShell\\;gsElasticity\\;${{ github.event.repository.name }}' -Q \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index cf3d4d0..19ec531 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,10 +46,10 @@ else() add_subdirectory(tutorials EXCLUDE_FROM_ALL) endif(GISMO_BUILD_EXAMPLES) -# # add unittests -# aux_gs_cpp_directory(${PROJECT_SOURCE_DIR}/unittests unittests_SRCS) -# set(gismo_UNITTESTS ${gismo_UNITTESTS} ${unittests_SRCS} -# CACHE INTERNAL "gismo list of unittests") +# add unittests +aux_gs_cpp_directory(${PROJECT_SOURCE_DIR}/unittests unittests_SRCS) +set(gismo_UNITTESTS ${gismo_UNITTESTS} ${unittests_SRCS} + CACHE INTERNAL "gismo list of unittests") # needed?: set(EXECUTABLE_OUTPUT_PATH ${CMAKE_BINARY_DIR}/bin/) diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt index 06dcbf9..8688172 100644 --- a/benchmarks/CMakeLists.txt +++ b/benchmarks/CMakeLists.txt @@ -11,7 +11,7 @@ aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR} FILES) foreach(file ${FILES}) add_gismo_executable(${file}) get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}-benchmarks") if(GISMO_BUILD_EXAMPLES) set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") else(GISMO_BUILD_EXAMPLES) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 689a545..6342ca0 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -11,7 +11,7 @@ aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR} FILES) foreach(file ${FILES}) add_gismo_executable(${file}) get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}-examples") if(GISMO_BUILD_EXAMPLES) set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") else(GISMO_BUILD_EXAMPLES) diff --git a/solvers/CMakeLists.txt b/solvers/CMakeLists.txt index 01fe917..13be6b8 100644 --- a/solvers/CMakeLists.txt +++ b/solvers/CMakeLists.txt @@ -11,7 +11,7 @@ aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR} FILES) foreach(file ${FILES}) add_gismo_executable(${file}) get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}-solvers") if(GISMO_BUILD_EXAMPLES) set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") else(GISMO_BUILD_EXAMPLES) diff --git a/src/gsStaticSolvers/gsStaticOpt.h b/src/gsStaticSolvers/gsStaticOpt.h index dd224a9..6322972 100644 --- a/src/gsStaticSolvers/gsStaticOpt.h +++ b/src/gsStaticSolvers/gsStaticOpt.h @@ -15,11 +15,27 @@ #include #pragma once - namespace gismo { -template +/** + * @file gsStaticOpt.h + * @brief Header file for the gsStaticOpt and gsOptProblemStatic classes. + * + * This file contains the definitions for the gsStaticOpt and gsOptProblemStatic classes, + * which are used for static optimization in structural analysis. + */ + +/** + * @class gsOptProblemStatic + * @brief A class representing a static optimization problem. + * + * @tparam T The coefficient type. + * + * This class inherits from gsOptProblem and provides functionality for evaluating + * the objective function and its gradient for static optimization problems. + */ +template class gsOptProblemStatic : public gsOptProblem { protected: @@ -29,7 +45,13 @@ class gsOptProblemStatic : public gsOptProblem typedef typename gsStructuralAnalysisOps::ALResidual_t ALResidual_t; public: - + /** + * @brief Constructor for gsOptProblemStatic. + * + * @param residualFun The residual function. + * @param L A coefficient (unused but needs to be assigned). + * @param numDesignVars The number of design variables. + */ gsOptProblemStatic(const typename gsStructuralAnalysisOps::Residual_t & residualFun, const T & L, index_t numDesignVars) : m_residualFun(residualFun), @@ -40,6 +62,13 @@ class gsOptProblemStatic : public gsOptProblem m_curDesign.setZero(); } + /** + * @brief Constructor for gsOptProblemStatic with arc-length residual function. + * + * @param ALresidualFun The arc-length residual function. + * @param L A coefficient. + * @param numDesignVars The number of design variables. + */ gsOptProblemStatic(const typename gsStructuralAnalysisOps::ALResidual_t & ALresidualFun, const T & L, index_t numDesignVars) : m_ALresidualFun(ALresidualFun), @@ -55,7 +84,12 @@ class gsOptProblemStatic : public gsOptProblem } public: - + /** + * @brief Evaluates the objective function. + * + * @param u The input vector. + * @return The norm of the residual. + */ T evalObj( const gsAsConstVector & u ) const { gsVector result; @@ -63,6 +97,12 @@ class gsOptProblemStatic : public gsOptProblem return result.norm(); } + /** + * @brief Computes the gradient of the objective function. + * + * @param u The input vector. + * @param result The output vector for the gradient. + */ void gradObj_into( const gsAsConstVector & u, gsAsVector & result) const { gsVector tmp; @@ -71,10 +111,6 @@ class gsOptProblemStatic : public gsOptProblem result = -tmp; } - // void evalCon_into( const gsAsConstVector & u, gsAsVector & result) const; - - // void jacobCon_into( const gsAsConstVector & u, gsAsVector & result) const; - private: Residual_t m_residualFun; @@ -98,10 +134,15 @@ class gsOptProblemStatic : public gsOptProblem }; /** - * @brief Static solver using the Dynamic Relaxation method + * @class gsStaticOpt + * @brief Static solver using the Dynamic Relaxation method. * - * @tparam T coefficient type + * @tparam T The coefficient type. + * @tparam Optimizer The optimizer type (default is gsGradientDescent). * + * This class inherits from gsStaticBase and provides functionality for solving + * static optimization problems using the Dynamic Relaxation method. + * \ingroup gsStaticBase */ template > @@ -117,9 +158,10 @@ class gsStaticOpt : public gsStaticBase public: /** - * @brief Constructor + * @brief Constructor for gsStaticOpt. * - * @param[in] Residual The residual + * @param Residual The residual function. + * @param numDofs The number of degrees of freedom. */ gsStaticOpt( const Residual_t &Residual, const index_t numDofs ) : @@ -131,9 +173,10 @@ class gsStaticOpt : public gsStaticBase } /** - * @brief Constructs a new instance. + * @brief Constructor for gsStaticOpt with arc-length residual function. * - * @param[in] ALResidual The residual as arc-length object + * @param ALResidual The arc-length residual function. + * @param numDofs The number of degrees of freedom. */ gsStaticOpt( const ALResidual_t & ALResidual, const index_t numDofs ) : @@ -146,6 +189,11 @@ class gsStaticOpt : public gsStaticBase } public: + /** + * @brief Returns the optimizer options. + * + * @return A reference to the optimizer options. + */ gsOptionList & optimizerOptions() { return m_optimizer.options(); } /// gsStaticBase base functions @@ -158,7 +206,7 @@ class gsStaticOpt : public gsStaticBase // /// See \ref gsStaticBase // void initOutput() override; - + // /// See \ref gsStaticBase // void stepOutput(index_t k) override; @@ -183,22 +231,12 @@ class gsStaticOpt : public gsStaticBase GISMO_NO_IMPLEMENTATION; } -// Class-specific functions -// public: -// void initialize(); - protected: /// See \ref solve() void _solve(); /// Initializes the method void _init(); -// gsVector _computeResidual(const gsVector & U); - -// public: - -// /// Return the residual norm -// T residualNorm() const { return m_R.norm(); } protected: Optimizer m_optimizer; diff --git a/src/gsStaticSolvers/gsStaticOpt.hpp b/src/gsStaticSolvers/gsStaticOpt.hpp index 7ca2c6a..101e0ea 100644 --- a/src/gsStaticSolvers/gsStaticOpt.hpp +++ b/src/gsStaticSolvers/gsStaticOpt.hpp @@ -93,150 +93,11 @@ gsStatus gsStaticOpt::solve() { this->getOptions(); m_optimizer.solve(m_U+m_DeltaU); - + m_DeltaU = m_optimizer.currentDesign() - m_U; m_U += m_DeltaU; m_numIterations = m_optimizer.iterations(); return gsStatus::Success; } - - - -//! [OptProblem] - -// template -// void gsStaticOpt::getOptions() -// { -// Base::getOptions(); -// } - -// template -// gsStatus gsStaticOpt::solve() -// { -// try -// { -// _solve(); -// m_status = gsStatus::Success; -// } -// catch (int errorCode) -// { -// if (errorCode==1) -// m_status = gsStatus::NotConverged; -// else if (errorCode==2) -// m_status = gsStatus::AssemblyError; -// else if (errorCode==3) -// m_status = gsStatus::SolverError; -// else -// m_status = gsStatus::OtherError; -// } -// catch (...) -// { -// m_status = gsStatus::OtherError; -// } -// return m_status; -// } - -// template -// void gsStaticOpt::_solve() -// { -// m_Eks.clear(); -// m_Eks.reserve(m_maxIterations); - -// if (m_verbose) initOutput(); - -// _start(); - -// m_Ek0 = m_Ek; -// m_Eks.push_back(m_Ek); - -// if (m_verbose != 0) stepOutput(0); -// index_t resetIt = 0; -// for (m_numIterations=1; m_numIterations!=m_maxIterations; m_numIterations++, resetIt++) -// { -// _iteration(); -// if ((m_c==0 && m_Ek_prev > m_Ek) || resetIt==m_resetIterations)// || (m_Ek/m_Ek_prev > 1/m_tolE && m_Ek_prev!=0)) -// { -// resetIt = 0; -// _peak(); -// } - -// if (m_verbose!=0) -// if (m_numIterations % m_verbose == 0 || m_verbose==-1 ) stepOutput(m_numIterations); - -// m_Eks.push_back(m_Ek); - -// m_residualOld = m_residual; - -// if (m_residual/m_residualIni < m_tolF && m_Ek/m_Ek0 < m_tolE && m_deltaU.norm()/m_DeltaU.norm() < m_tolU) -// { -// m_U += m_DeltaU; -// gsDebug <<"Converged: \n"; -// gsDebug <<"\t |R|/|R0| = "< -// void gsStaticOpt::initialize() -// { -// this->reset(); -// getOptions(); -// } - -// template -// void gsStaticOpt::reset() -// { -// Base::reset(); -// m_status = gsStatus::NotStarted; -// } - -// template -// void gsStaticOpt::_init() -// { -// this->reset(); -// defaultOptions(); -// } - -// template -// gsVector gsStaticOpt::_computeResidual(const gsVector & U) -// { -// gsVector resVec; -// if (!m_residualFun(U, resVec)) -// throw 2; -// return resVec; -// } - -// template -// gsOptProblemStatic::gsOptProblemStatic() -// { - -// // Number of design variables -// m_numDesignVars = 2; -// // Number of constraints -// m_numConstraints = 0; -// // Number of non-zeros in the Jacobian of the constraints -// m_numConJacNonZero = 0; - -// m_desLowerBounds.resize(m_numDesignVars); -// m_desUpperBounds.resize(m_numDesignVars); - -// m_desLowerBounds.setConstant(-1e19); -// m_desUpperBounds.setConstant( 1e19); - -// // we initialize x in bounds, in the upper right quadrant -// m_curDesign.resize(m_numDesignVars,1); -// m_curDesign.setZero(); -// } - } // namespace gismo diff --git a/tutorials/CMakeLists.txt b/tutorials/CMakeLists.txt index 54267b8..a5d6550 100644 --- a/tutorials/CMakeLists.txt +++ b/tutorials/CMakeLists.txt @@ -11,7 +11,7 @@ aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR} FILES) foreach(file ${FILES}) add_gismo_executable(${file}) get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}-tutorials") if(GISMO_BUILD_EXAMPLES) set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") else(GISMO_BUILD_EXAMPLES) diff --git a/unittests/gsStaticSolver_test.cpp b/unittests/gsStaticSolver_test.cpp new file mode 100644 index 0000000..6cb0378 --- /dev/null +++ b/unittests/gsStaticSolver_test.cpp @@ -0,0 +1,421 @@ +/** @file gsThinShellAssembler_test.cpp + + @brief Provides unittests for the gsThinShellAssembler class + + * Balloon: unit-test based on a hyperelastic balloon inflated by a follower pressure. + This test allows to test the follower pressure as well as the hyperelastic material models (incompressible) + + * UAT: unit-test based on a uni-axial tension test + This test allows to test the hyperelastic material models (incompressible and compressible) + + * Modal: unit-test based on a modal analysis + This test allows to test the mass matrix and the compressible material model + + + == BASIC REFERENCE == + - TEST(NAME_OF_TEST) { body_of_test } + - TEST_FIXTURE(NAME_OF_FIXTURE,NAME_OF_TEST){ body_of_test } + + == CHECK MACRO REFERENCE == + - CHECK(EXPR); + - CHECK_EQUAL(EXPECTED,ACTUAL); + - CHECK_CLOSE(EXPECTED,ACTUAL,EPSILON); + - CHECK_ARRAY_EQUAL(EXPECTED,ACTUAL,LENGTH); + - CHECK_ARRAY_CLOSE(EXPECTED,ACTUAL,LENGTH,EPSILON); + - CHECK_ARRAY2D_EQUAL(EXPECTED,ACTUAL,ROWCOUNT,COLCOUNT); + - CHECK_ARRAY2D_CLOSE(EXPECTED,ACTUAL,ROWCOUNT,COLCOUNT,EPSILON); + - CHECK_THROW(EXPR,EXCEPTION_TYPE_EXPECTED); + + == TIME CONSTRAINTS == + - UNITTEST_TIME_CONSTRAINT(TIME_IN_MILLISECONDS); + - UNITTEST_TIME_CONSTRAINT_EXEMPT(); + + == MORE INFO == + See: https://unittest-cpp.github.io/ + + Author(s): H.M.Verhelst (2019 - ..., TU Delft) + **/ + +#include "gismo_unittest.h" // Brings in G+Smo and the UnitTest++ framework +#ifdef gsKLShell_ENABLED +#include +#endif + +#ifdef gsOptim_ENABLED +#include +#endif + +#include +#include +#include +#include + +SUITE(gsThinShellAssembler_test) // The suite should have the same name as the file +{ + +#ifdef gsKLShell_ENABLED + // solver: 0: newton, 1: DR, 2: OPT + std::pair UAT_analytical(const std::vector solver, const index_t material=1, const index_t impl=1, const bool Compressibility=false); + std::pair UAT_numerical(const std::vector solver, const index_t material=1, const index_t impl=1, const bool Compressibility=false); + void UAT_CHECK(const std::vector solver, const index_t material=1, const index_t impl=1, const bool Compressibility=false); +#endif + +///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +#ifdef gsKLShell_ENABLED + TEST(StaticSolver_UAT_NR) + { + std::vector solver({0}); + UAT_CHECK(solver); + } + + // TEST(StaticSolver_UAT_DR) + // { + // std::vector solver({1}); + // UAT_CHECK(solver); + // } + +#ifdef gsOptim_ENABLED + TEST(StaticSolver_UAT_OP) + { + std::vector solver({2}); + UAT_CHECK(solver); + } +#endif + +#ifdef gsOptim_ENABLED + TEST(StaticSolver_UAT_OP_NR) + { + std::vector solver({2,0}); + UAT_CHECK(solver); + } +#endif + + TEST(StaticSolver_UAT_DR_NR) + { + std::vector solver({1,0}); + UAT_CHECK(solver); + } + +#ifdef gsOptim_ENABLED + TEST(StaticSolver_UAT_DR_OP) + { + std::vector solver({1,2}); + UAT_CHECK(solver); + } +#endif + + std::pair UAT_numerical(const std::vector solver, const index_t material, const index_t impl, const bool Compressibility) + { + //! [Parse command line] + index_t numRefine = 1; + index_t numElevate = 1; + + real_t E_modulus = 1.0; + real_t PoissonRatio; + real_t Density = 1.0; + real_t Ratio = 7.0; + + real_t mu = 1.5e6; + real_t thickness = 0.001; + + real_t alpha1,alpha2,alpha3,mu1,mu2,mu3; + alpha1 = 1.3; + mu1 = 6.3e5/4.225e5*mu; + alpha2 = 5.0; + mu2 = 0.012e5/4.225e5*mu; + alpha3 = -2.0; + mu3 = -0.1e5/4.225e5*mu; + + if (!Compressibility) + PoissonRatio = 0.5; + else + PoissonRatio = 0.45; + + E_modulus = 2*mu*(1+PoissonRatio); + + //! [Parse command line] + + //! [Read input file] + gsMultiPatch<> mp, mp_def; + + mp.addPatch( gsNurbsCreator<>::BSplineSquare(1) ); // degree + + if (numElevate!=0) + mp.degreeElevate(numElevate); + + // h-refine + for (int r =0; r < numRefine; ++r) + mp.uniformRefine(); + + mp_def = mp; + + //! [Refinement] + gsMultiBasis<> dbasis(mp); + + gsBoundaryConditions<> bc; + bc.setGeoMap(mp); + + gsPointLoads pLoads = gsPointLoads(); + + real_t lambda = 2.0; + gsConstantFunction<> displx(lambda-1.0,2); + + GISMO_ENSURE(mp.targetDim()==2,"Geometry must be planar (targetDim=2)!"); + bc.addCondition(boundary::west, condition_type::dirichlet, 0, 0, false, 0 ); + + bc.addCondition(boundary::east, condition_type::dirichlet, &displx, 0, false, 0 ); + + bc.addCondition(boundary::south, condition_type::dirichlet, 0, 0, false, 1 ); + + //! [Refinement] + + // Linear isotropic material model + gsVector<> tmp(2); + tmp.setZero(); + gsConstantFunction<> force(tmp,2); + gsFunctionExpr<> t(std::to_string(thickness),2); + gsFunctionExpr<> E(std::to_string(E_modulus),2); + gsFunctionExpr<> nu(std::to_string(PoissonRatio),2); + gsFunctionExpr<> rho(std::to_string(Density),2); + gsConstantFunction<> ratio(Ratio,2); + + gsConstantFunction<> alpha1fun(alpha1,2); + gsConstantFunction<> mu1fun(mu1,2); + gsConstantFunction<> alpha2fun(alpha2,2); + gsConstantFunction<> mu2fun(mu2,2); + gsConstantFunction<> alpha3fun(alpha3,2); + gsConstantFunction<> mu3fun(mu3,2); + + std::vector*> parameters(3); + parameters[0] = &E; + parameters[1] = ν + parameters[2] = ∶ + gsMaterialMatrixBase* materialMatrix; + + if (material==4) + { + parameters.resize(8); + parameters[0] = &E; + parameters[1] = ν + parameters[2] = &mu1fun; + parameters[3] = &alpha1fun; + parameters[4] = &mu2fun; + parameters[5] = &alpha2fun; + parameters[6] = &mu3fun; + parameters[7] = &alpha3fun; + } + + gsOptionList options; + if (material==0) + { + GISMO_ERROR("This test is not available for SvK models"); + } + else + { + options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",material); + options.addSwitch("Compressibility","Compressibility: (false): Imcompressible | (true): Compressible",Compressibility); + options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",impl); + materialMatrix = getMaterialMatrix<2,real_t>(mp,t,parameters,rho,options); + } + + gsThinShellAssemblerBase* assembler; + assembler = new gsThinShellAssembler<2, real_t, false >(mp,dbasis,bc,force,materialMatrix); + + assembler->setPointLoads(pLoads); + + // Function for the Jacobian + typedef std::function (gsVector const &)> Jacobian_t; + typedef std::function (gsVector const &) > Residual_t; + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&mp_def](gsVector const &x, gsSparseMatrix & m) + { + ThinShellAssemblerStatus status; + assembler->constructSolution(x,mp_def); + status = assembler->assembleMatrix(mp_def); + m = assembler->matrix(); + return status == ThinShellAssemblerStatus::Success; + }; + + gsStructuralAnalysisOps::Residual_t Residual = [&assembler,&mp_def](gsVector const &x, gsVector & result) + { + ThinShellAssemblerStatus status; + assembler->constructSolution(x,mp_def); + status = assembler->assembleVector(mp_def); + result = assembler->rhs(); + return status == ThinShellAssemblerStatus::Success; + }; + + // Define Matrices + assembler->assemble(); + gsSparseMatrix<> K = assembler->matrix(); + gsVector<> F = assembler->rhs(); + assembler->assembleMass(true); + gsVector<> M = assembler->rhs(); + + gsStaticDR DRM(M,F,Residual); + gsOptionList DROptions = DRM.options(); + DROptions.setReal("damping",0.0); + DROptions.setReal("alpha",1e9); + DROptions.setInt("maxIt",1e6); + DROptions.setReal("tolF",1e-8); + DROptions.setReal("tolU",1e-6); + DROptions.setReal("tolE",1e-4); + DROptions.setInt("verbose",0); + DRM.setOptions(DROptions); + DRM.initialize(); + +#ifdef gsOptim_ENABLED + gsStaticOpt::LBFGS> OPT(Residual,assembler->numDofs()); +#else + gsStaticOpt> OPT(Residual,assembler->numDofs()); +#endif + gsOptionList OPTOptions = OPT.options(); + OPTOptions.setInt("maxIt",100); + OPTOptions.setReal("tolF",1e-8); + OPTOptions.setReal("tolU",1e-6); + OPTOptions.setInt("verbose",0); + OPT.setOptions(OPTOptions); + OPT.initialize(); + + gsStaticNewton NWT(K,F,Jacobian,Residual); + gsOptionList NWTOptions = NWT.options(); + NWTOptions.setInt("maxIt",20); + NWTOptions.setReal("tolF",1e-8); + NWTOptions.setReal("tolU",1e-6); + NWTOptions.setInt("verbose",0); + NWT.setOptions(NWTOptions); + NWT.initialize(); + + std::vector*> solverVec(solver.size()); + for (size_t s=0; s!=solver.size(); s++) + if (solver[s]==0) + solverVec[s] = &NWT; + else if (solver[s]==1) + solverVec[s] = &DRM; + else if (solver[s]==2) + solverVec[s] = &OPT; + + gsStaticComposite compositeSolver(solverVec); + compositeSolver.initialize(); + compositeSolver.solve(); + CHECK(compositeSolver.status() == gsStatus::Success); + mp_def = assembler->constructSolution(compositeSolver.solution()); + + gsMultiPatch<> deformation = mp_def; + for (size_t k = 0; k != mp_def.nPatches(); ++k) + deformation.patch(k).coefs() -= mp.patch(k).coefs(); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // Check solutions + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // NOTE: all analytical solutions for compressible materials are fixed for displ=1; (lambda=2) + + // Compute stretches (should be the same everywhere) + // Ordering: lambda(0) < lambda(1); lambda(2) is ALWAYS the through-thickness stretch + gsVector<> pt(2); + pt<<1,0; + gsMatrix<> lambdas = assembler->computePrincipalStretches(pt,mp_def,0); + + // Get the total force on the tension boundary + patchSide ps(0,boundary::east); + gsMatrix<> forceVector = assembler->boundaryForce(mp_def,ps); + real_t sideForce = forceVector.sum(); + real_t S = -sideForce / (thickness*lambdas(0)*lambdas(2)); + real_t L = lambdas(0); + + std::pair result; + result.first = L; + result.second = S; + + delete materialMatrix; + delete assembler; + + return result; + } + + std::pair UAT_analytical(const std::vector solver, const index_t material, const index_t impl, const bool Compressibility) + { + real_t PoissonRatio; + real_t Ratio = 7.0; + + real_t mu = 1.5e6; + + real_t alpha1,alpha2,alpha3,mu1,mu2,mu3; + alpha1 = 1.3; + mu1 = 6.3e5/4.225e5*mu; + alpha2 = 5.0; + mu2 = 0.012e5/4.225e5*mu; + alpha3 = -2.0; + mu3 = -0.1e5/4.225e5*mu; + + if (!Compressibility) + PoissonRatio = 0.5; + else + PoissonRatio = 0.45; + + real_t lambda = 2.0; + + real_t San,J,K,Lan; + if (material==1 && Compressibility) + { + K = 2*mu*(1+PoissonRatio)/(3-6*PoissonRatio); + J = 1.105598565;// specific for lambda==2!! + San = lambda*(0.5*mu*(-(2*(math::pow(lambda,2)+2*J/lambda))/(3*math::pow(J,2./3.)*lambda)+2*lambda/math::pow(J,2./3.))+0.25*K*(2*math::pow(J,2)/lambda-2./lambda))/J; + Lan = math::pow(J/lambda,0.5); + } + else if (material==1 && !Compressibility) + { + San = mu * (lambda*lambda - 1/lambda); + Lan = math::pow(1./lambda,0.5); + } + else if (material==3 && Compressibility) + { + real_t c2 = 1.0 / (Ratio+1); + real_t c1 = 1.0 - c2; + K = 2*mu*(1+PoissonRatio)/(3-6*PoissonRatio); + J = 1.099905842;// specific for lambda==2!! + San = lambda*(0.5*c1*mu*(-(2*(math::pow(lambda,2)+2*J/lambda))/(3*math::pow(J,2./3.)*lambda)+2*lambda/math::pow(J,2./3.))+0.5*c2*mu*(-(4*(2*lambda*J+math::pow(J,2)/math::pow(lambda,2)))/(3*math::pow(J,4./3.)*lambda)+4/math::pow(J,1./3.))+0.25*K*(2*math::pow(J,2)/lambda-2/lambda))/J; + Lan = math::pow(J/lambda,0.5); + } + else if (material==3 && !Compressibility) + { + real_t c2 = 1.0 / (Ratio+1); + real_t c1 = 1.0 - c2; + San =-mu*(c2*lambda*lambda+c2/lambda+c1)/lambda+lambda*(c1*lambda*mu+2*c2*mu); + Lan = math::pow(1./lambda,0.5); + } + else if (material==4 && Compressibility) + { + K = 2*mu*(1+PoissonRatio)/(3-6*PoissonRatio); + J = 1.088778638;// specific for lambda==2!! + San = 1./J* (lambda *( mu1*(2*math::pow(lambda/math::pow(J,1./3.),alpha1)*alpha1/(3*lambda)-2*math::pow(math::pow(J/lambda,0.5)/math::pow(J,1./3.),alpha1)*alpha1/(3*lambda))/alpha1+mu2*(2*math::pow(lambda/math::pow(J,1./3.),alpha2)*alpha2/(3*lambda)-2*math::pow(math::pow(J/lambda,0.5)/math::pow(J,1./3.),alpha2)*alpha2/(3*lambda))/alpha2+mu3*(2*math::pow(lambda/math::pow(J,1./3.),alpha3)*alpha3/(3*lambda)-2*math::pow(math::pow(J/lambda,0.5)/math::pow(J,1./3.),alpha3)*alpha3/(3*lambda))/alpha3+0.25*K*(2*math::pow(J,2)/lambda-2/lambda) ) ); + Lan = math::pow(J/lambda,0.5); + } + else if (material==4 && !Compressibility) + { + San =-mu1*math::pow((1./lambda),0.5*alpha1)-mu2*math::pow((1./lambda),0.5*alpha2)-mu3*math::pow((1./lambda),0.5*alpha3)+mu1*math::pow(lambda,alpha1)+mu2*math::pow(lambda,alpha2)+mu3*math::pow(lambda,alpha3); + Lan = math::pow(1./lambda,0.5); + } + else + GISMO_ERROR("Material not treated"); + + std::pair result; + result.first = Lan; + result.second = San; + return result; + } + + void UAT_CHECK(const std::vector solver, const index_t material, const index_t impl, const bool Compressibility) + { + if (material==4 && impl!=3) + CHECK(true); + + real_t Lnum, Snum, Lana, Sana; + std::tie(Lnum,Snum) = UAT_numerical(solver,material,impl,Compressibility); + std::tie(Lana,Sana) = UAT_analytical(solver,material,impl,Compressibility); + CHECK_CLOSE(std::abs(Lnum-Lana)/Lana,0,1e-7); + } +#endif + +} From 3233d571ca3cd603e261c445dab0658206b0ec27 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Thu, 2 Jan 2025 20:29:59 +0100 Subject: [PATCH 06/29] add loop for failure of load step --- benchmarks/benchmark_cylinder_DC.cpp | 29 +++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 2fb4d10..fe116c0 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -63,7 +63,7 @@ int main (int argc, char** argv) index_t maxitOPT = 1000; bool OP = false; // * Newton Raphson - index_t maxitNR = 20; + index_t maxitNR = 50; bool NR = false; @@ -110,7 +110,6 @@ int main (int argc, char** argv) gsOptionList opts; fd.getFirst(opts); - real_t R = 0.25, H = 1.0; real_t thickness = 0.05e-3; real_t YoungsModulus = 1e9; real_t PoissonRatio = 0.5; @@ -341,14 +340,32 @@ int main (int argc, char** argv) file<updateBCs(BCs); StaticSolver.setDisplacement(U); StaticSolver.solve(); + + if (StaticSolver.status() != gsStatus::Success) + { + dL = dL/2; + D = D - dL; + continue; + } + else + { + dL = std::min(dL0,1-D); + D = D + dL; + k++; + } + U = StaticSolver.solution(); // 1. Get all the coefficients (including the ones from the eliminated BCs.) @@ -550,12 +567,6 @@ int main (int argc, char** argv) << D << "," << Moment << "\n"; file.close(); - - - - - // controlDC.step(dL); - D+= dL; } if (plot) From bed82fe89156d3eabf3e1ade1f18cabd19fcd57b Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Fri, 3 Jan 2025 09:59:38 +0100 Subject: [PATCH 07/29] incremental in z --- benchmarks/benchmark_cylinder_DC.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index fe116c0..f0549aa 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -159,8 +159,7 @@ int main (int argc, char** argv) dz_val = 125e-3; else if (testCase==1) dz_val = 1.0; - - gsConstantFunction<> dz(dz_val,3); + gsFunctionExpr<> dz(util::to_string(dz_val) + "*t",3); // Boundary conditions gsBoundaryConditions<> BCs; @@ -349,6 +348,7 @@ int main (int argc, char** argv) gsInfo<<"Displacement step "<updateBCs(BCs); StaticSolver.setDisplacement(U); StaticSolver.solve(); From 8b8a050a49e3958225ee3e0347b81fb53d8ce55a Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 15 Jan 2025 08:53:38 +0100 Subject: [PATCH 08/29] update benchmark --- benchmarks/benchmark_cylinder_DC.cpp | 33 ++++++++++++++-------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index f0549aa..4f95ab6 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -55,7 +55,7 @@ int main (int argc, char** argv) real_t dL = 1; // Solvers // * Dynamic Relaxation - index_t maxitDR = 1e4; + index_t maxitDR = 1e5; real_t alpha = 1.0; real_t damping = 0.1; bool DR = false; @@ -159,7 +159,7 @@ int main (int argc, char** argv) dz_val = 125e-3; else if (testCase==1) dz_val = 1.0; - gsFunctionExpr<> dz(util::to_string(dz_val) + "*t",3); + gsFunctionExpr<> dz(util::to_string(dz_val) + "*1",3); // Boundary conditions gsBoundaryConditions<> BCs; @@ -203,8 +203,8 @@ int main (int argc, char** argv) parameters[0] = &E; parameters[1] = ν - gsMaterialMatrixBase * materialMatrix; - gsMaterialMatrixBase * materialMatrixTFT; + gsMaterialMatrixBase::uPtr materialMatrix; + gsMaterialMatrixBase::uPtr materialMatrixTFT; gsThinShellAssemblerBase* assembler; gsOptionList options; @@ -215,16 +215,17 @@ int main (int argc, char** argv) gsMaterialMatrixContainer materialMatrixContainer; if (TFT) { - materialMatrixTFT = new gsMaterialMatrixTFT<3,real_t,true>(materialMatrix); + materialMatrixTFT = memory::make_unique(new gsMaterialMatrixTFT<3,real_t,true>(*materialMatrix)); + materialMatrixTFT->options().setReal("SlackMultiplier",1e-6); for (size_t p = 0; p!=mp.nPatches(); p++) - materialMatrixContainer.add(materialMatrixTFT); + materialMatrixContainer.add(*materialMatrixTFT); // materialMatrixTFT->options().setReal("SlackMultiplier",1e-6); assembler = new gsThinShellAssembler<3, real_t, false >(geom,dbasis,BCs,force,materialMatrixContainer); } else { for (size_t p = 0; p!=mp.nPatches(); p++) - materialMatrixContainer.add(materialMatrix); + materialMatrixContainer.add(*materialMatrix); assembler = new gsThinShellAssembler<3, real_t, true >(geom,dbasis,BCs,force,materialMatrixContainer); } @@ -284,8 +285,9 @@ int main (int argc, char** argv) DynamicRelaxationOptions.setReal("damping",damping); DynamicRelaxationOptions.setReal("alpha",alpha); DynamicRelaxationOptions.setInt("maxIt",maxitDR); - DynamicRelaxationOptions.setReal("tolE",1e-4); - DynamicRelaxationOptions.setReal("tol",1e-2); + DynamicRelaxationOptions.setReal("tolE",1e-2); + DynamicRelaxationOptions.setReal("tolF",1e-4); + DynamicRelaxationOptions.setReal("tolU",1e-4); DynamicRelaxationOptions.setInt("verbose",1); DynamicRelaxationSolver.setOptions(DynamicRelaxationOptions); DynamicRelaxationSolver.initialize(); @@ -293,8 +295,8 @@ int main (int argc, char** argv) gsStaticOpt::LBFGS> OptimizerSolver(Residual,assembler->numDofs()); gsOptionList OptimizerSolverOptions = OptimizerSolver.options(); OptimizerSolverOptions.setInt("maxIt",maxitOPT); - OptimizerSolverOptions.setReal("tolF",1e-4); - OptimizerSolverOptions.setReal("tolU",1e-4); + OptimizerSolverOptions.setReal("tolF",1e-5); + OptimizerSolverOptions.setReal("tolU",1e-5); OptimizerSolverOptions.setInt("verbose",1); OptimizerSolver.setOptions(OptimizerSolverOptions); OptimizerSolver.initialize(); @@ -304,8 +306,8 @@ int main (int argc, char** argv) gsOptionList NewtonOptions = NewtonSolver.options(); NewtonOptions.setInt("verbose",true); NewtonOptions.setInt("maxIt",maxitNR); - NewtonOptions.setReal("tolU",1e-4); - NewtonOptions.setReal("tolF",1e-6); + NewtonOptions.setReal("tolU",1e-6); + NewtonOptions.setReal("tolF",1e-8); NewtonSolver.setOptions(NewtonOptions); std::vector *> solvers; @@ -575,10 +577,7 @@ int main (int argc, char** argv) collectionTF.save(); } - delete materialMatrix; delete assembler; - if (TFT) - delete materialMatrixTFT; return EXIT_SUCCESS; } @@ -590,4 +589,4 @@ int main(int argc, char *argv[]) gsWarn<<"G+Smo is not compiled with the gsKLShell module."; return EXIT_FAILURE; } -#endif \ No newline at end of file +#endif From 4f749d5fafa4481bc31b95d398cb95f1b47262cc Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 15 Jan 2025 11:51:45 +0100 Subject: [PATCH 09/29] remove slack multiplier --- benchmarks/benchmark_cylinder_DC.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 4f95ab6..5950db4 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -216,10 +216,8 @@ int main (int argc, char** argv) if (TFT) { materialMatrixTFT = memory::make_unique(new gsMaterialMatrixTFT<3,real_t,true>(*materialMatrix)); - materialMatrixTFT->options().setReal("SlackMultiplier",1e-6); for (size_t p = 0; p!=mp.nPatches(); p++) materialMatrixContainer.add(*materialMatrixTFT); - // materialMatrixTFT->options().setReal("SlackMultiplier",1e-6); assembler = new gsThinShellAssembler<3, real_t, false >(geom,dbasis,BCs,force,materialMatrixContainer); } else From ad3bb576759e1dd5a6770e1bb732ebfb348354ab Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 15 Jan 2025 12:18:29 +0100 Subject: [PATCH 10/29] change output moment --- solvers/static_shell_multipatch_XML_torsion.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/solvers/static_shell_multipatch_XML_torsion.cpp b/solvers/static_shell_multipatch_XML_torsion.cpp index c639613..b433840 100644 --- a/solvers/static_shell_multipatch_XML_torsion.cpp +++ b/solvers/static_shell_multipatch_XML_torsion.cpp @@ -560,10 +560,10 @@ int main (int argc, char** argv) auto dEf= -( deriv2(v,sn(deff).normalized().tr() ) + deriv2(deff,var1(v,deff) ) ) * reshape(m2,3,3); - gsDebugVar(evA.integral(( + gsInfo<<"Moment = "<<(evA.integral(( N*dEm.tr() + M*dEf.tr() - )*meas(ori))); + )*meas(ori)))<<"\n"; // std::vector bContainer({patchSide(0,boundary::north),patchSide(1,boundary::north),patchSide(2,boundary::north),patchSide(3,boundary::north)}); @@ -714,4 +714,4 @@ int main(int argc, char *argv[]) gsWarn<<"G+Smo is not compiled with the gsKLShell module."; return EXIT_FAILURE; } -#endif \ No newline at end of file +#endif From dea276042b12bf3e3056f8462d40df9f0b0071a0 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 15 Jan 2025 15:45:09 +0100 Subject: [PATCH 11/29] small cleanup --- benchmarks/benchmark_cylinder_DC.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 5950db4..83d1a7e 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -218,7 +218,7 @@ int main (int argc, char** argv) materialMatrixTFT = memory::make_unique(new gsMaterialMatrixTFT<3,real_t,true>(*materialMatrix)); for (size_t p = 0; p!=mp.nPatches(); p++) materialMatrixContainer.add(*materialMatrixTFT); - assembler = new gsThinShellAssembler<3, real_t, false >(geom,dbasis,BCs,force,materialMatrixContainer); + assembler = new gsThinShellAssembler<3, real_t, false>(geom,dbasis,BCs,force,materialMatrixContainer); } else { @@ -256,7 +256,7 @@ int main (int argc, char** argv) }; // Function for the Residual - gsStructuralAnalysisOps::Residual_t Residual = [&dx,&dy,&BCs,&assembler,&bb2, &geom](gsVector const &x, gsVector & result) + gsStructuralAnalysisOps::Residual_t Residual = [&assembler,&bb2, &geom](gsVector const &x, gsVector & result) { ThinShellAssemblerStatus status; gsMatrix solFull = assembler->fullSolutionVector(x); From 0e70bca93a61810fd8f8d84ac2e42982d3ab98c1 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 15 Jan 2025 15:56:35 +0100 Subject: [PATCH 12/29] staticDR will re-scale residual correctly in second iteration, iff residual==0 in first iteration --- src/gsStaticSolvers/gsStaticDR.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/gsStaticSolvers/gsStaticDR.hpp b/src/gsStaticSolvers/gsStaticDR.hpp index 10240a9..f904038 100644 --- a/src/gsStaticSolvers/gsStaticDR.hpp +++ b/src/gsStaticSolvers/gsStaticDR.hpp @@ -188,6 +188,7 @@ void gsStaticDR::_iteration() m_Ek_prev = m_Ek; m_R = _computeResidual(m_U+m_DeltaU) - m_damp.cwiseProduct(m_v); m_residual = m_R.norm(); + if (m_residualIni==0) m_residualIni = m_residual; //---------------------------------------------------------------------------------- m_v += m_dt * m_massInv.cwiseProduct(m_R); // Velocities at t+dt/2 m_deltaU = m_dt * m_v; // Velocities at t+dt @@ -231,7 +232,7 @@ void gsStaticDR::_start() // m_residual = m_R.norm(); m_residual = m_forcing.norm(); // If the residual is 0 (e.g. with purely displacment loading), we set it to 1 to allow divisions - if (m_residual==0) m_residual=1; + if (m_residual==0) m_residual=0; // All residual norms are equal m_residualIni = m_residualOld = m_residual; } From c568540eafc59dd16f467ea2e013132dcb302c26 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Mon, 27 Jan 2025 21:16:38 +0100 Subject: [PATCH 13/29] add skip after three refinements --- benchmarks/benchmark_cylinder_DC.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 83d1a7e..3473b09 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -326,6 +326,8 @@ int main (int argc, char** argv) } gsStaticComposite StaticSolver(solvers); + StaticSolver.options().setInt("verbose",1); + StaticSolver.initialize(); // gsControlDisplacement controlDC(&StaticSolver); // Displacement-controlled simulation real_t dL0 = dL; @@ -353,7 +355,7 @@ int main (int argc, char** argv) StaticSolver.setDisplacement(U); StaticSolver.solve(); - if (StaticSolver.status() != gsStatus::Success) + if (StaticSolver.status() != gsStatus::Success && dL/dL0 > math::pow(0.5,3)) { dL = dL/2; D = D - dL; From 7dadeda846fbf859addd17caa01b5715ac339000 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Mon, 27 Jan 2025 22:10:49 +0100 Subject: [PATCH 14/29] fix stepping --- benchmarks/benchmark_cylinder_DC.cpp | 29 +++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 3473b09..496e41b 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -52,6 +52,7 @@ int main (int argc, char** argv) int step = 10; bool TFT = false; + bool membrane= false; real_t dL = 1; // Solvers // * Dynamic Relaxation @@ -84,6 +85,7 @@ int main (int argc, char** argv) // Material settings cmd.addSwitch("TFT", "Use TFT material matrix", TFT); + cmd.addSwitch("membrane", "Use membrane material", membrane); // Test case cmd.addInt("t", "test", "Test case: (0): Annulus | (1): Cylinder", testCase); @@ -220,6 +222,12 @@ int main (int argc, char** argv) materialMatrixContainer.add(*materialMatrixTFT); assembler = new gsThinShellAssembler<3, real_t, false>(geom,dbasis,BCs,force,materialMatrixContainer); } + else if (membrane) + { + for (size_t p = 0; p!=mp.nPatches(); p++) + materialMatrixContainer.add(*materialMatrix); + assembler = new gsThinShellAssembler<3, real_t, false >(geom,dbasis,BCs,force,materialMatrixContainer); + } else { for (size_t p = 0; p!=mp.nPatches(); p++) @@ -342,11 +350,8 @@ int main (int argc, char** argv) file<<0<<","<<0<<"\n"; file.close(); index_t k=0; - while (D <= 1) + while (true) { - if (gsClose(D,1.0,1e-8)) - break; - gsInfo<<"Displacement step "< math::pow(0.5,3)) + if ((StaticSolver.status() != gsStatus::Success && dL/dL0 > math::pow(0.5,3)) || gsClose(D,0.4,1e-8)) { dL = dL/2; D = D - dL; continue; } - else - { - dL = std::min(dL0,1-D); - D = D + dL; - k++; - } U = StaticSolver.solution(); @@ -569,6 +568,14 @@ int main (int argc, char** argv) << D << "," << Moment << "\n"; file.close(); + + if (gsClose(D,1.0,1e-8)) + break; + + // GO TO THE NEXT STEPs + dL = std::min(dL0,1-D); + D = D + dL; + k++; } if (plot) From c68ffb4bd9cd260ab26c65d219cee18d98c2dace Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Tue, 28 Jan 2025 09:47:16 +0100 Subject: [PATCH 15/29] small --- benchmarks/benchmark_cylinder_DC.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 496e41b..0cd047d 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -182,6 +182,8 @@ int main (int argc, char** argv) GISMO_ERROR("Test case "< Date: Tue, 28 Jan 2025 09:53:39 +0100 Subject: [PATCH 16/29] add output line, remove intermediate artificial step --- benchmarks/benchmark_cylinder_DC.cpp | 39 +++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 0cd047d..4fcd620 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -40,6 +40,27 @@ #include using namespace gismo; +void writeToCSVfile(std::string name, gsMatrix<> matrix) +{ + std::ofstream file(name.c_str()); + for(int i = 0; i < matrix.rows(); i++) + { + for(int j = 0; j < matrix.cols(); j++) + { + std::string str = std::to_string(matrix(i,j)); + if(j+1 == matrix.cols()) + { + file< math::pow(0.5,3)) || gsClose(D,0.4,1e-8)) + if ((StaticSolver.status() != gsStatus::Success && dL/dL0 > math::pow(0.5,3))) { dL = dL/2; D = D - dL; @@ -405,6 +426,22 @@ int main (int argc, char** argv) fileNameTF = gsFileManager::getFilename(fileNameTF); for (size_t p=0; p coords(2,N), supp, result; + index_t dir = 1; + for (index_t patch=0; patch!=def.nPieces(); patch++) + { + supp = def.piece(patch).support(); + coords.row(dir).setConstant( 0.5*( supp(dir,1)-supp(dir,0) )+supp(dir,0) ); + coords.row(1-dir).setLinSpaced(N,0,1); + coords.row(1-dir).array() *= ( supp(1-dir,1)-supp(1-dir,0) ); + coords.row(1-dir).array() += supp(1-dir,0); + def.piece(patch).eval_into(coords,result); + result.transposeInPlace(); + writeToCSVfile(dirname + "/" + "line_step=" + std::to_string(k) + "_dir=" + std::to_string(dir) + ".csv",result); + } } ///////////////////////////////////////////////////////////////////////////////// From 0299b859aadfae61bd60e576304996ea17eac514 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Thu, 30 Jan 2025 14:05:23 +0100 Subject: [PATCH 17/29] add lines per coordinate --- benchmarks/benchmark_cylinder_DC.cpp | 30 ++++++++++++++++++---------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 4fcd620..92ad337 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -431,18 +431,26 @@ int main (int argc, char** argv) index_t N = 100; gsMatrix<> coords(2,N), supp, result; index_t dir = 1; - for (index_t patch=0; patch!=def.nPieces(); patch++) - { - supp = def.piece(patch).support(); - coords.row(dir).setConstant( 0.5*( supp(dir,1)-supp(dir,0) )+supp(dir,0) ); - coords.row(1-dir).setLinSpaced(N,0,1); - coords.row(1-dir).array() *= ( supp(1-dir,1)-supp(1-dir,0) ); - coords.row(1-dir).array() += supp(1-dir,0); - def.piece(patch).eval_into(coords,result); - result.transposeInPlace(); - writeToCSVfile(dirname + "/" + "line_step=" + std::to_string(k) + "_dir=" + std::to_string(dir) + ".csv",result); + std::vector par; + if (testCase==0) + par = {0.5,0.7,0.9}; + else if (testCase==1) + par = {0.5,0.75}; + else + GISMO_ERROR("Test case "< Date: Fri, 31 Jan 2025 17:25:14 +0100 Subject: [PATCH 18/29] improved export of lines --- benchmarks/benchmark_cylinder_DC.cpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 92ad337..ae4fc45 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -429,7 +429,7 @@ int main (int argc, char** argv) index_t N = 100; - gsMatrix<> coords(2,N), supp, result; + gsMatrix<> coords(2,N), supp, result, tmp; index_t dir = 1; std::vector par; if (testCase==0) @@ -438,7 +438,10 @@ int main (int argc, char** argv) par = {0.5,0.75}; else GISMO_ERROR("Test case "< Date: Sat, 1 Feb 2025 15:48:26 +0100 Subject: [PATCH 19/29] add checkpoints --- benchmarks/benchmark_cylinder_DC.cpp | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index ae4fc45..b255dcc 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -373,6 +373,12 @@ int main (int argc, char** argv) file<<0<<","<<0<<"\n"; file.close(); index_t k=0; + std::queue Dcheckpoints; + Dcheckpoints.push(0.7); + Dcheckpoints.push(0.8); + Dcheckpoints.push(0.9); + Dcheckpoints.push(1.0); + dL = std::min(dL,Dcheckpoints.front()); while (true) { gsInfo<<"Displacement step "< Date: Sat, 1 Feb 2025 17:12:40 +0100 Subject: [PATCH 20/29] add membrane option to Pillow --- benchmarks/benchmark_Pillow.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/benchmarks/benchmark_Pillow.cpp b/benchmarks/benchmark_Pillow.cpp index 5c36db3..b6a9bf4 100644 --- a/benchmarks/benchmark_Pillow.cpp +++ b/benchmarks/benchmark_Pillow.cpp @@ -42,6 +42,7 @@ int main (int argc, char** argv) bool DR = false; bool NR = false; int verbose = 0; + bool membrane = false; index_t Compressibility = 1; index_t material = 1; @@ -86,6 +87,7 @@ int main (int argc, char** argv) cmd.addSwitch("split", "Split to a multi-patch", split); cmd.addSwitch("smooth", "Use a smooth basis (maximum regularity)", smooth); cmd.addInt("v","verbose", "0: no; 1: iteration output; 2: Full matrix and vector output", verbose); + cmd.addSwitch("membrane", "Membrane element", membrane); try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } GISMO_ASSERT(NR || DR,"At least a Newton Raphson or a Dynamic Relaxation solver needs to be enabled. Run with option NR or DR."); @@ -236,7 +238,11 @@ int main (int argc, char** argv) // materialMatrixTFT->options().setSwitch("Explicit",true); gsThinShellAssemblerBase* assembler; - if (TFT) + if (!membrane && TFT) + assembler = new gsThinShellAssembler<3, real_t, true >(mp,dbasis,BCs,force,materialMatrixTFT); + else if (!membrane && !TFT) + assembler = new gsThinShellAssembler<3, real_t, true >(mp,dbasis,BCs,force,materialMatrix); + else if (membrane && TFT) assembler = new gsThinShellAssembler<3, real_t, false >(mp,dbasis,BCs,force,materialMatrixTFT); else assembler = new gsThinShellAssembler<3, real_t, false >(mp,dbasis,BCs,force,materialMatrix); From 09ed5a22c5f371286b58396ed0bdd27a8ee326ba Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 25 Jun 2025 11:08:42 +0200 Subject: [PATCH 21/29] add tolerance for DynamicRelaxation --- benchmarks/benchmark_Pillow.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/benchmarks/benchmark_Pillow.cpp b/benchmarks/benchmark_Pillow.cpp index b6a9bf4..f116d49 100644 --- a/benchmarks/benchmark_Pillow.cpp +++ b/benchmarks/benchmark_Pillow.cpp @@ -44,6 +44,8 @@ int main (int argc, char** argv) int verbose = 0; bool membrane = false; + real tolDR = 1e-1; + index_t Compressibility = 1; index_t material = 1; index_t impl = 1; // 1= analytical, 2= generalized, 3= spectral @@ -79,6 +81,8 @@ int main (int argc, char** argv) cmd.addReal( "L", "maxLoad", "Maximum load", maxLoad ); + cmd.addReal( "T", "tolDR", "Tolerance for the DR solver", tolDR ); + cmd.addSwitch("plot", "Plot result in ParaView format", plot); cmd.addSwitch("write", "Write data to file", write); cmd.addSwitch("stress", "Plot stress in ParaView format", stress); @@ -311,7 +315,7 @@ int main (int argc, char** argv) DROptions.setReal("damping",damping); DROptions.setReal("alpha",alpha); DROptions.setInt("maxIt",1e6); - DROptions.setReal("tol",1e-2); + DROptions.setReal("tol",tolDR); DROptions.setReal("tolE",1e-5); DROptions.setInt("verbose",verbose); DROptions.setInt("ResetIt",(index_t)(100)); @@ -443,6 +447,8 @@ int main (int argc, char** argv) gsVector<> data(1); data< Date: Wed, 25 Jun 2025 11:27:07 +0200 Subject: [PATCH 22/29] add comments to torsion example --- benchmarks/benchmark_cylinder_DC.cpp | 242 ++++++++++----------------- 1 file changed, 88 insertions(+), 154 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index b255dcc..4492a42 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -359,16 +359,13 @@ int main (int argc, char** argv) gsStaticComposite StaticSolver(solvers); StaticSolver.options().setInt("verbose",1); StaticSolver.initialize(); - // gsControlDisplacement controlDC(&StaticSolver); - // Displacement-controlled simulation - real_t dL0 = dL; - real_t D = dL; + gsMatrix<> U = gsMatrix<>::Zero(assembler->numDofs(),1); std::string writeName = dirname + "/out.csv"; std::ofstream file; file.open(writeName); - file<<"Load factor, Moment\n"; + file<<"Load factor, Moment1, Moment2, Moment3\n"; file< A(1,1); - // A.setIntegrationElements(dbasis); - // gsExprEvaluator<> evA(A); - // gsExprAssembler<>::space u = A.getSpace(bb2,3); - - // // gsConstantFunction<> one(1.0,3); - // // gsBoundaryConditions<> bc_dummy; - // // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - - // gsFunctionExpr<> Mz_x("sqrt(x^2+y^2)*cos(atan2(y,x))",3); - // gsFunctionExpr<> Mz_y("sqrt(x^2+y^2)*sin(atan2(y,x))",3); - // gsBoundaryConditions<> bc_dummy; - // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - - // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - - - // bc_dummy.setGeoMap(geom); - // u.setup(bc_dummy,dirichlet::interpolation,-1); - // gsDebugVar(u.mapper()); - - // gsMatrix<> coefs(bb2.size(),3); - // for (size_t d=0; d!=geom.geoDim(); d++) - // for (index_t k=0; k!=bb2.size(); k++) - // { - // const index_t ii = u.mapper().index(k,0,d); - // if (u.mapper().is_free(k,0,d)) - // coefs(k,d) = 0.0; - // else if (u.mapper().is_boundary(k,0,d)) - // coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); - // else - // GISMO_ERROR("WHAT?"); - // } - - // gsMappedSpline<2,real_t> geo(bb2,coefs); - - - gsExprAssembler<> A(1,1); - A.setIntegrationElements(dbasis); - gsExprEvaluator<> evA(A); - gsExprAssembler<>::space u = A.getSpace(dbasis,3); - - // gsConstantFunction<> one(1.0,3); // for X reaction force - // gsBoundaryConditions<> bc_dummy; - // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - - gsBoundaryConditions<> bc_dummy; - gsFunctionExpr<> Mz_x("sqrt(x^2+y^2)*cos(atan2(y,x))",3); - gsFunctionExpr<> Mz_y("-sqrt(x^2+y^2)*sin(atan2(y,x))",3); // X loads contribute negatively to M_z, so we need to flip the sign - bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - - bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - - /* - gsBoundaryConditions<> bc_dummy; - gsMultiPatch<> mp_x = mp.coord(0); - gsMultiPatch<> mp_y = mp.coord(1); - // X loads contribute negatively to M_z, so we need to flip the sign - for (size_t p=0; p!=mp_y.nPatches(); p++) - mp_y.patch(p).coefs().array() *= -1; - bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); - bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); - bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); - bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); - - bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); - bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); - bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); - bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); - */ - bc_dummy.setGeoMap(mp); - u.setup(bc_dummy,dirichlet::interpolation,-1); - gsDebugVar(u.mapper()); - - gsMultiPatch<> geo; - for (size_t b=0; b!=dbasis.nBases(); b++) + // Define vectors representing the distance from the origin to the boundary for the moment around the x and y axes. + std::vector> Mz_x_vec(3); + std::vector> Mz_y_vec(3); + Mz_x_vec[0] = gsFunctionExpr<>("sqrt(x^2+y^2)*cos(atan2(y,x))",3); + Mz_x_vec[1] = dx; + Mz_x_vec[2] = gsFunctionExpr<>("sqrt(x^2+y^2)*cos(atan2(y,x)+pi/2*t)",3); + Mz_x_vec[2].set_t(D); + + Mz_y_vec[0] = gsFunctionExpr<>("-sqrt(x^2+y^2)*sin(atan2(y,x))",3); + Mz_y_vec[1] = dy; + Mz_y_vec[2] = gsFunctionExpr<>("sqrt(x^2+y^2)*sin(atan2(y,x)+pi/2*t)",3); + Mz_y_vec[2].set_t(D); + + std::vector Moments(3); + + for (index_t k=0; k!=Mz_x_vec.size(); k++) { - gsMatrix<> coefs(dbasis.basis(b).size(),mp.geoDim()); - coefs.setZero(); - for (size_t d=0; d!=mp.geoDim(); d++) - for (index_t k=0; k!=dbasis.basis(b).size(); k++) + gsExprAssembler<> A(1,1); + A.setIntegrationElements(dbasis); + gsExprEvaluator<> evA(A); + gsExprAssembler<>::space u = A.getSpace(dbasis,3); + + // Define dummy boundary conditions. + // These conditions are used to define the test space on the boundary and compute its values there. + gsBoundaryConditions<> bc_dummy; + bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_y_vec[k], 0 ,false,0); + bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_y_vec[k], 0 ,false,0); + bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_y_vec[k], 0 ,false,0); + bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_y_vec[k], 0 ,false,0); + + bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_x_vec[k], 0 ,false,1); + bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_x_vec[k], 0 ,false,1); + bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_x_vec[k], 0 ,false,1); + bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_x_vec[k], 0 ,false,1); + + bc_dummy.setGeoMap(mp); + u.setup(bc_dummy,dirichlet::interpolation,-1); + + // Set all coefficients corresponding to the space to zero, except for the boundary conditions. + gsMultiPatch<> geo; + for (size_t b=0; b!=dbasis.nBases(); b++) { - const index_t ii = u.mapper().index(k,b,d); - if (u.mapper().is_free(k,b,d)) - coefs(k,d) = 0.0; - else if (u.mapper().is_boundary(k,b,d)) - coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); - else - GISMO_ERROR("WHAT?"); + gsMatrix<> coefs(dbasis.basis(b).size(),mp.geoDim()); + coefs.setZero(); + for (size_t d=0; d!=mp.geoDim(); d++) + for (index_t k=0; k!=dbasis.basis(b).size(); k++) + { + const index_t ii = u.mapper().index(k,b,d); + if (u.mapper().is_free(k,b,d)) + coefs(k,d) = 0.0; + else if (u.mapper().is_boundary(k,b,d)) + coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); + else + GISMO_ERROR("WHAT?"); + } + geo.addPatch(*dbasis.basis(b).makeGeometry(give(coefs))); } - geo.addPatch(*dbasis.basis(b).makeGeometry(give(coefs))); + // Cleanup the assembler + A.initSystem(); + // Set the geometry map and the deformation map. + gsExprAssembler<>::geometryMap ori = A.getMap(mp); + gsExprAssembler<>::geometryMap deff = A.getMap(def); + // Register the custom solution with coefficients defined above to the assembler. + auto v = A.getCoeff(geo); + // Add the stress tensors based on the materials defined for the shell solver + gsMaterialMatrixIntegrate S0f(materialMatrixContainer,&mp,&def); + gsMaterialMatrixIntegrate S1f(materialMatrixContainer,&mp,&def); + auto S0 = A.getCoeff(S0f); + auto S1 = A.getCoeff(S1f); + // Helper matrix for flexural components (see gsThinShellAssembler::assembleMatrix) + gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2); + auto m2 = A.getCoeff(mult2t); + // Define the components used in the variational formulation, using the solution v as the deformation + auto N = S0.tr(); + auto dEm= flat( jac(deff).tr() * jac(v) ) ; + auto M = S1.tr(); + auto dEf= -( deriv2(v,sn(deff).normalized().tr() ) + deriv2(deff,var1(v,deff) ) ) * reshape(m2,3,3); + // Integrate the variational formulation, this gives the moments. + Moments[k] = evA.integral(( + N*dEm.tr() + + M*dEf.tr() + )*meas(ori) + ); } - // for (size_t b=0; b!=dbasis.nBases(); b++) - // for (size_t d=0; d!=2; d++) - // for (index_t k=0; k!=dbasis.basis(b).size(); k++) - // { - // const index_t ii = u.mapper().index(k,0,d); - // if (u.mapper().is_free(k,0,d)) - // coefs(k,d) = 0.0; - // else if (u.mapper().is_boundary(k,0,d)) - // coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); - // else - // GISMO_ERROR("WHAT?"); - // } - - // gsGeometry<>::uPtr geo = dbasis.basis(0).makeGeometry(give(coefs)); - - A.initSystem(); - - gsExprAssembler<>::geometryMap ori = A.getMap(mp); - gsExprAssembler<>::geometryMap deff = A.getMap(def); - - auto v = A.getCoeff(geo); - - gsMaterialMatrixIntegrate S0f(materialMatrixContainer,&mp,&def); - gsMaterialMatrixIntegrate S1f(materialMatrixContainer,&mp,&def); - auto S0 = A.getCoeff(S0f); - auto S1 = A.getCoeff(S1f); - - // Helper matrix for flexural components - gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2); - auto m2 = A.getCoeff(mult2t); - - auto N = S0.tr(); - auto dEm= flat( jac(deff).tr() * jac(v) ) ; - auto M = S1.tr(); - auto dEf= -( deriv2(v,sn(deff).normalized().tr() ) + deriv2(deff,var1(v,deff) ) ) * reshape(m2,3,3); - - real_t Moment = evA.integral(( - N*dEm.tr() - + M*dEf.tr() - )*meas(ori) - ); - file.open(writeName,std::ofstream::out | std::ofstream::app); file<< std::setprecision(10) << D << "," - << Moment << "\n"; + << Moments[0] << "," + << Moments[1] << "," + << Moments[2] << "\n"; file.close(); if (gsClose(D,Dcheckpoints.front(),1e-8)) From 10554c4a988fe1d79461983f01ba2862d5bf75e6 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 25 Jun 2025 13:44:48 +0200 Subject: [PATCH 23/29] add lumping --- solvers/static_shell_XML.cpp | 9 +- .../static_shell_multipatch_XML_torsion.cpp | 717 ------------------ 2 files changed, 4 insertions(+), 722 deletions(-) delete mode 100644 solvers/static_shell_multipatch_XML_torsion.cpp diff --git a/solvers/static_shell_XML.cpp b/solvers/static_shell_XML.cpp index 44fcb0d..7857898 100644 --- a/solvers/static_shell_XML.cpp +++ b/solvers/static_shell_XML.cpp @@ -240,9 +240,8 @@ int main (int argc, char** argv) assembler->assemble(); gsSparseMatrix<> K = assembler->matrix(); gsVector<> F = assembler->rhs(); - assembler->assembleMass(false); - gsSparseMatrix<> Mmat = assembler->matrix(); - gsVector<> M = Mmat.toDense().rowwise().sum(); + assembler->assembleMass(true); + gsVector<> M = assembler->rhs(); // Function for the Jacobian gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&mp_def](gsVector const &x, gsSparseMatrix & m) @@ -330,7 +329,7 @@ int main (int argc, char** argv) result.col(k) = deformation.patch(refPatches.at(k)).eval(refPoints.col(k)); writer.add(result); - gsWrite(mp_def,"deformed"); + gsWrite(mp_def,dirname + sep + "deformed"); } // ! [Export visualization in ParaView] @@ -392,7 +391,7 @@ int main (int argc, char** argv) gsWriteParaview(stretchDir1,dirname + sep + "PrincipalDirection1",5000); gsWriteParaview(stretchDir2,dirname + sep + "PrincipalDirection2",5000); gsWriteParaview(stretchDir3,dirname + sep + "PrincipalDirection3",5000); - + } delete assembler; diff --git a/solvers/static_shell_multipatch_XML_torsion.cpp b/solvers/static_shell_multipatch_XML_torsion.cpp deleted file mode 100644 index b433840..0000000 --- a/solvers/static_shell_multipatch_XML_torsion.cpp +++ /dev/null @@ -1,717 +0,0 @@ -/** @file static_shell_multipatch_XML.cpp - - @brief Blackbox solver for static shell analysis using unstructured splines - - 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 (2019-..., TU Delft) -*/ - -#include - -#ifdef gsKLShell_ENABLED -#include -#include -#include -#endif - -#include -#include -#include - -#include - -#ifdef gsUnstructuredSplines_ENABLED -#include -#include -#include -#endif - -#include - -using namespace gismo; - -#ifdef gsKLShell_ENABLED -#ifdef gsUnstructuredSplines_ENABLED -int main (int argc, char** argv) -{ - // Input options - int numElevate = 0; - int numRefine = 1; - bool plot = false; - bool write = false; - bool stress = false; - bool DR = false; - bool NR = false; - - bool membrane = false; - - index_t method = 0; - - real_t perturb = 0; - - // Arc length method options - - std::string bvp; - - std::string dirname = "."; - - std::string wn("data.csv"); - - gsCmdLine cmd("Shell static solver for multipatches."); - - cmd.addInt("r","hRefine", "Number of dyadic h-refinement (bisection) steps to perform before solving", numRefine); - cmd.addInt("e","degreeElevation", "Number of degree elevation steps to perform on the Geometry's basis before solving", numElevate); - cmd.addInt("m","method", "Smoothing method to use: 0: smoothInterfaces, 1: Almost C1, 2: D-Patch", method); - - cmd.addReal( "P" , "perturb" , "Perturb the z coordinate of the refined geometry", perturb ); - - cmd.addSwitch("membrane", "Membrane element", membrane); - - cmd.addSwitch("plot", "Plot result in ParaView format", plot); - cmd.addSwitch("write", "Write data to file", write); - cmd.addSwitch("stress", "Plot stress in ParaView format", stress); - cmd.addSwitch("DR", "Use Dynamic Relaxation", DR); - cmd.addSwitch("NR", "Use Newton Raphson", NR); - - cmd.addString("i","inputFile", "Input file", bvp); - cmd.addString("o","outputDir", "Output directory", dirname); - - - try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } - - // ![Initialize members] - gsMultiPatch<> mp; - // Boundary conditions - gsBoundaryConditions<> BCs; - - gsFunctionExpr<> forceFun; - gsFunctionExpr<> pressFun; - - gsPointLoads pLoads = gsPointLoads(); - gsMatrix<> points,loads; - gsMatrix pid_ploads; - - gsOptionList DROptions, NROptions, assemblerOptions; - // ![Initialize members] - - // ![Read data] - gsFileData fd(bvp); - if (fd.hasAny>()) - { - gsInfo<<"Reading geometry from "<(path + geomFileName,mp); - } - gsInfo<<"Finished\n"; - - // p-refine - for (size_t p=0; p!=mp.nPatches(); p++) - { - for(index_t i = 0; i< numElevate; ++i) - { - if (dynamic_cast * >(&mp.patch(p))) - { - gsWarn<<"Degree elevation applied"<<"\n"; - mp.patch(p).degreeElevate(); // Elevate the degree - } - else - mp.patch(p).degreeIncrease(); // Elevate the degree - } - - // h-refine - for(index_t i = 0; i< numRefine; ++i) - mp.patch(p).uniformRefine(); - } - - if (perturb!=0) - for (size_t p=0; p!=mp.nPatches(); p++) - { - mp.patch(p).coefs().col(2).setRandom(); - mp.patch(p).coefs().col(2) *= perturb; - } - - gsMultiBasis<> dbasis(mp,true); - gsInfo<<"Basis (patch 0): "<< mp.patch(0).basis() << "\n"; - if (gsTensorNurbsBasis<2,real_t> * basis = dynamic_cast * >(&mp.basis(0))) - for (index_t dim = 0; dim!=2; dim++) - gsDebug<<"dir "<knots(dim).asMatrix()<<"\n"; - - gsInfo<<"Looking for material matrices ...\n"; - gsMaterialMatrixContainer materialMatrixContainer; - gsMaterialMatrixBase * materialMatrix; - if (fd.hasAny>()) - { - gsInfo<<"Reading material matrix container (ID=11) ..."; - fd.getId(11,materialMatrixContainer); - } - else - { - gsInfo<<"Reading material matrix (ID=10) ..."; - materialMatrix = fd.getId>(10).release(); - for (size_t p = 0; p!=mp.nPatches(); p++) - materialMatrixContainer.add(materialMatrix); - } - gsInfo<<"Finished\n"; - gsInfo<<"Finished\n"; - - gsInfo<<"Reading boundary conditions (ID=20) ..."; - fd.getId(20,BCs); - gsInfo<<"Finished\n"; - - gsInfo<<"Reading force function (ID=21) ..."; - fd.getId(21,forceFun); - gsInfo<<"Finished\n"; - bool pressure = false; - if ( fd.hasId(22) ) - { - gsInfo<<"Reading pressure function (ID=22) ..."; - pressure = true; - fd.getId(22,pressFun); - gsInfo<<"Finished\n"; - } - // Point loads - gsInfo<<"Reading point load locations (ID=30) ..."; - if ( fd.hasId(30) ) fd.getId(30,points); - gsInfo<<"Finished\n"; - gsInfo<<"Reading point load vectors (ID=31) ..."; - if ( fd.hasId(31) ) fd.getId(31,loads); - gsInfo<<"Finished\n"; - gsInfo<<"Reading point load patch indices (ID=32) ..."; - if ( fd.hasId(32) ) fd.getId(32,pid_ploads); - gsInfo<<"Finished\n"; - - if ( !fd.hasId(30) || !fd.hasId(31) || !fd.hasId(32) ) - pid_ploads = gsMatrix::Zero(1,points.cols()); - - for (index_t k =0; k!=points.cols(); k++) - pLoads.addLoad(points.col(k), loads.col(k), pid_ploads.at(k) ); // in parametric domain! - - - // Reference points - gsMatrix refPatches; - gsMatrix<> refPoints, refPars, refValue; // todo: add refValue.. - gsInfo<<"Reading reference point locations (ID=50) ..."; - if ( fd.hasId(50) ) fd.getId(50,refPoints); - gsInfo<<"Finished\n"; - gsInfo<<"Reading reference patches (ID=51) ..."; - if ( fd.hasId(51) ) fd.getId(51,refPatches); - gsInfo<<"Finished\n"; - gsInfo<<"Reading reference values (ID=52) ..."; - if ( fd.hasId(52) ) fd.getId(52,refValue); - gsInfo<<"Finished\n"; - - if ( !fd.hasId(50) || !fd.hasId(51) || !fd.hasId(52) ) - refValue = gsMatrix<>::Zero(mp.geoDim(),refPoints.cols()); - GISMO_ENSURE(refPatches.cols()==refPoints.cols(),"Number of reference points and patches do not match"); - - if (refPoints.rows()==2) - { - refPars = refPoints; - gsInfo<<"Reference points are provided in parametric coordinates.\n"; - } - else if (refPoints.rows()==3) - gsInfo<<"Reference points are provided in physical coordinates.\n"; - else - gsInfo<<"No reference points are provided.\n"; - - gsInfo<<"Reading DR solver options (ID=90) ..."; - if ( fd.hasId(90) ) fd.getId(90,DROptions); - gsInfo<<"Finished\n"; - gsInfo<<"Reading NR solver options (ID=91) ..."; - if ( fd.hasId(91) ) fd.getId(91,NROptions); - gsInfo<<"Finished\n"; - gsInfo<<"Reading assembler options (ID=92) ..."; - if ( fd.hasId(92) ) fd.getId(92,assemblerOptions); - gsInfo<<"Finished\n"; - - // ![Read data] - - // Fix path - dirname = gsFileManager::getCanonicRepresentation(dirname,true); - char sep = gsFileManager::getNativePathSeparator(); - gsFileManager::mkdir(dirname); - - // plot geometry - if (plot) - gsWriteParaview(mp,"mp",1000,true); - - mp.computeTopology(); - - // Make unstructured spline - gsMultiPatch<> geom; - gsMappedBasis<2,real_t> bb2; - gsSparseMatrix<> global2local; - if (method==0) - { - gsSmoothInterfaces<2,real_t> smoothInterfaces(mp); - smoothInterfaces.options().setSwitch("SharpCorners",false); - smoothInterfaces.compute(); - smoothInterfaces.matrix_into(global2local); - - global2local = global2local.transpose(); - geom = smoothInterfaces.exportToPatches(); - dbasis = smoothInterfaces.localBasis(); - } - else if (method==1) - { - gsAlmostC1<2,real_t> almostC1(mp); - almostC1.options().setSwitch("SharpCorners",true); - almostC1.compute(); - almostC1.matrix_into(global2local); - - global2local = global2local.transpose(); - geom = almostC1.exportToPatches(); - dbasis = almostC1.localBasis(); - } - else if (method==2) - { - gsDPatch<2,real_t> dpatch(mp); - dpatch.options().setSwitch("SharpCorners",true); - dpatch.compute(); - dpatch.matrix_into(global2local); - - global2local = global2local.transpose(); - geom = dpatch.exportToPatches(); - dbasis = dpatch.localBasis(); - } - else - GISMO_ERROR("Method "<* assembler; - if (membrane) - assembler = new gsThinShellAssembler<3, real_t, false >(geom,dbasis,BCs,forceFun,materialMatrixContainer); - else - assembler = new gsThinShellAssembler<3, real_t, true >(geom,dbasis,BCs,forceFun,materialMatrixContainer); - assembler->setOptions(assemblerOptions); - assembler->options().setInt("Continuity",-1); - assembler->setSpaceBasis(bb2); - assembler->setPointLoads(pLoads); - if (pressure) - assembler->setPressure(pressFun); - - // Assemble linear system to obtain the force vector - gsDebugVar(assembler->numDofs()); - assembler->assemble(); - gsSparseMatrix<> K = assembler->matrix(); - gsVector<> F = assembler->rhs(); - assembler->assembleMass(true); - gsVector<> M = assembler->rhs(); - - gsStructuralAnalysisOps::Jacobian_t Jacobian; - gsStructuralAnalysisOps::Residual_t Residual; - // Function for the Jacobian - Jacobian = [&assembler,&bb2,&geom](gsVector const &x, gsSparseMatrix & m) - { - ThinShellAssemblerStatus status; - gsMatrix solFull = assembler->fullSolutionVector(x); - size_t d = geom.targetDim(); - GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); - solFull.resize(solFull.rows()/d,d); - gsMappedSpline<2,real_t> mspline(bb2,solFull); - gsFunctionSum def(&geom,&mspline); - assembler->assembleMatrix(def); - status = assembler->assembleMatrix(def); - m = assembler->matrix(); - return status == ThinShellAssemblerStatus::Success; - }; - // Function for the Residual - Residual = [&assembler,&bb2,&geom](gsVector const &x, gsVector & result) - { - ThinShellAssemblerStatus status; - gsMatrix solFull = assembler->fullSolutionVector(x); - size_t d = geom.targetDim(); - GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); - solFull.resize(solFull.rows()/d,d); - - gsMappedSpline<2,real_t> mspline(bb2,solFull); - gsFunctionSum def(&geom,&mspline); - - status = assembler->assembleVector(def); - result = assembler->rhs(); // assembler rhs - force = Finternal - return status == ThinShellAssemblerStatus::Success; - }; - - gsStructuralAnalysisOutput writer(dirname + sep + wn,refPoints); - std::vector pointheaders = {"x","y","z"}; - std::vector otherheaders = {}; - writer.init(pointheaders,otherheaders); - - gsStaticNewton LIN(K,F); - LIN.setOptions(NROptions); - - gsStaticDR DRM(M,F,Residual); - DRM.setOptions(DROptions); - - gsStaticNewton NRM(K,F,Jacobian,Residual); - NRM.setOptions(NROptions); - - std::vector *> solvers{&LIN}; - if (DR) - { - gsInfo<<"Using Dynamic Relaxation solver\n"; - solvers.push_back(&DRM); - } - if (NR) - { - gsInfo<<"Using Newton-Raphson solver\n"; - solvers.push_back(&NRM); - } - gsStaticComposite solver(solvers); - solver.initialize(); - solver.solve(); - GISMO_ASSERT(solver.converged(),"Solver failed"); - gsVector<> solVector = solver.solution(); - - // ! [Export visualization in ParaView] - if (plot) - { - gsInfo<<"Plotting in Paraview...\n"; - std::string fileName; - /// Make a gsMappedSpline to represent the solution - // 1. Get all the coefficients (including the ones from the eliminated BCs.) - gsMatrix solFull = assembler->fullSolutionVector(solVector); - - // 2. Reshape all the coefficients to a Nx3 matrix - size_t d = geom.targetDim(); - GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); - solFull.resize(solFull.rows()/d,d); - - // 3. Make the mapped spline - gsMappedSpline<2,real_t> mspline(bb2,solFull); - - // 4. Create deformation spline - gsFunctionSum def(&geom,&mspline); - - // 5. Plot the mapped spline on the original geometry - gsPiecewiseFunction<> displacements; - assembler->constructStress(def,displacements,stress_type::displacement); - fileName = dirname + sep + "solution"; - gsWriteParaview(def,displacements,fileName,1000,"_"); - - // 5. Construct stress - gsPiecewiseFunction<> tensionFields; - assembler->constructStress(def,tensionFields,stress_type::tension_field); - fileName = dirname + sep + "tensionfield"; - gsWriteParaview(def,tensionFields,fileName,1000,"_"); - - - // gsExprAssembler<> A(1,1); - // A.setIntegrationElements(dbasis); - // gsExprEvaluator<> evA(A); - // gsExprAssembler<>::space u = A.getSpace(bb2,3); - - // // gsConstantFunction<> one(1.0,3); - // // gsBoundaryConditions<> bc_dummy; - // // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - - // gsFunctionExpr<> Mz_x("sqrt(x^2+y^2)*cos(atan2(y,x))",3); - // gsFunctionExpr<> Mz_y("sqrt(x^2+y^2)*sin(atan2(y,x))",3); - // gsBoundaryConditions<> bc_dummy; - // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - - // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - - - // bc_dummy.setGeoMap(geom); - // u.setup(bc_dummy,dirichlet::interpolation,-1); - // gsDebugVar(u.mapper()); - - // gsMatrix<> coefs(bb2.size(),3); - // for (size_t d=0; d!=geom.geoDim(); d++) - // for (index_t k=0; k!=bb2.size(); k++) - // { - // const index_t ii = u.mapper().index(k,0,d); - // if (u.mapper().is_free(k,0,d)) - // coefs(k,d) = 0.0; - // else if (u.mapper().is_boundary(k,0,d)) - // coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); - // else - // GISMO_ERROR("WHAT?"); - // } - - // gsMappedSpline<2,real_t> geo(bb2,coefs); - - - gsExprAssembler<> A(1,1); - A.setIntegrationElements(dbasis); - gsExprEvaluator<> evA(A); - gsExprAssembler<>::space u = A.getSpace(dbasis,3); - - // gsConstantFunction<> one(1.0,3); // for X reaction force - // gsBoundaryConditions<> bc_dummy; - // bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - // bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &one, 0 ,false,0); - - gsBoundaryConditions<> bc_dummy; - gsFunctionExpr<> Mz_x("sqrt(x^2+y^2)*cos(atan2(y,x))",3); - gsFunctionExpr<> Mz_y("-sqrt(x^2+y^2)*sin(atan2(y,x))",3); // X loads contribute negatively to M_z, so we need to flip the sign - bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_y, 0 ,false,0); - - bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, &Mz_x, 0 ,false,1); - -/* - gsBoundaryConditions<> bc_dummy; - gsMultiPatch<> mp_x = mp.coord(0); - gsMultiPatch<> mp_y = mp.coord(1); - // X loads contribute negatively to M_z, so we need to flip the sign - for (size_t p=0; p!=mp_y.nPatches(); p++) - mp_y.patch(p).coefs().array() *= -1; - bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); - bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); - bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); - bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, mp_y, 0 ,true,0); - - bc_dummy.addCondition(0,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); - bc_dummy.addCondition(1,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); - bc_dummy.addCondition(2,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); - bc_dummy.addCondition(3,boundary::north, condition_type::dirichlet, mp_x, 0 ,true,1); -*/ - bc_dummy.setGeoMap(mp); - u.setup(bc_dummy,dirichlet::interpolation,-1); - gsDebugVar(u.mapper()); - - gsMultiPatch<> geo; - for (size_t b=0; b!=dbasis.nBases(); b++) - { - gsMatrix<> coefs(dbasis.basis(b).size(),mp.geoDim()); - coefs.setZero(); - for (size_t d=0; d!=mp.geoDim(); d++) - for (index_t k=0; k!=dbasis.basis(b).size(); k++) - { - const index_t ii = u.mapper().index(k,b,d); - if (u.mapper().is_free(k,b,d)) - coefs(k,d) = 0.0; - else if (u.mapper().is_boundary(k,b,d)) - coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); - else - GISMO_ERROR("WHAT?"); - } - geo.addPatch(*dbasis.basis(b).makeGeometry(give(coefs))); - - } - - // for (size_t b=0; b!=dbasis.nBases(); b++) - // for (size_t d=0; d!=2; d++) - // for (index_t k=0; k!=dbasis.basis(b).size(); k++) - // { - // const index_t ii = u.mapper().index(k,0,d); - // if (u.mapper().is_free(k,0,d)) - // coefs(k,d) = 0.0; - // else if (u.mapper().is_boundary(k,0,d)) - // coefs(k,d) = u.fixedPart().at(u.mapper().global_to_bindex(ii)); - // else - // GISMO_ERROR("WHAT?"); - // } - - // gsGeometry<>::uPtr geo = dbasis.basis(0).makeGeometry(give(coefs)); - - A.initSystem(); - - gsExprAssembler<>::geometryMap ori = A.getMap(mp); - gsExprAssembler<>::geometryMap deff = A.getMap(def); - - auto v = A.getCoeff(geo); - - gsMaterialMatrixIntegrate S0f(materialMatrixContainer,&mp,&def); - gsMaterialMatrixIntegrate S1f(materialMatrixContainer,&mp,&def); - auto S0 = A.getCoeff(S0f); - auto S1 = A.getCoeff(S1f); - - // Helper matrix for flexural components - gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2); - auto m2 = A.getCoeff(mult2t); - - auto N = S0.tr(); - auto dEm= flat( jac(deff).tr() * jac(v) ) ; - auto M = S1.tr(); - auto dEf= -( deriv2(v,sn(deff).normalized().tr() ) + deriv2(deff,var1(v,deff) ) ) * reshape(m2,3,3); - - - gsInfo<<"Moment = "<<(evA.integral(( - N*dEm.tr() - + M*dEf.tr() - )*meas(ori)))<<"\n"; - // std::vector bContainer({patchSide(0,boundary::north),patchSide(1,boundary::north),patchSide(2,boundary::north),patchSide(3,boundary::north)}); - - - - - - - } - if (write) - { - // 1. Get all the coefficients (including the ones from the eliminated BCs.) - gsMatrix solFull = assembler->fullSolutionVector(solVector); - - // 2. Reshape all the coefficients to a Nx3 matrix - size_t d = geom.targetDim(); - GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); - solFull.resize(solFull.rows()/d,d); - - // 3. Make the mapped spline - gsMappedSpline<2,real_t> mspline(bb2,solFull); - - // 4. Make geom_def by projecting coefficients - gsDofMapper mapper(dbasis); - mapper.finalize(); - gsMatrix<> coefs; - gsMultiPatch<> geom_def = geom; - gsInfo<<"L2-Projection error of mspline on dbasis = "<::projectFunction(dbasis,mspline,geom,coefs)<<"\n"; - coefs.resize(coefs.rows()/geom.geoDim(),geom.geoDim()); - - index_t offset = 0; - for (size_t p = 0; p != geom_def.nPatches(); p++) - { - gsMatrix<> tmp_coefs = geom.patch(p).coefs(); - tmp_coefs += coefs.block(offset,0,mapper.patchSize(p),geom.geoDim()); - geom_def.patch(p) = give(*dbasis.basis(p).makeGeometry((tmp_coefs))); - offset += mapper.patchSize(p); - } - - // 5. Write the result - gsWrite(geom_def,dirname + sep + "deformed"); - - // gsMatrix<> result(mp.geoDim(),refPoints.cols()); - // for (index_t k=0; k!=refPoints.cols(); k++) - // result.col(k) = deformation.patch(refPatches.at(k)).eval(refPoints.col(k)); - // writer.add(result); - } - - // ! [Export visualization in ParaView] - if (stress) - { - gsField<> tensionField; - gsPiecewiseFunction<> tensionFields; - std::string fileName = dirname + sep + "tensionfield"; - - /// Make a gsMappedSpline to represent the solution - // 1. Get all the coefficients (including the ones from the eliminated BCs.) - gsMatrix solFull = assembler->fullSolutionVector(solVector); - - // 2. Reshape all the coefficients to a Nx3 matrix - size_t d = geom.targetDim(); - GISMO_ASSERT(solFull.rows() % d==0,"Rows of the solution vector does not match the number of control points"); - solFull.resize(solFull.rows()/d,d); - - // 3. Make the mapped spline - gsMappedSpline<2,real_t> mspline(bb2,solFull); - - // 4. Create deformation spline - gsFunctionSum def(&geom,&mspline); - - // 5. Construct stress - assembler->constructStress(def,tensionFields,stress_type::tension_field); - gsWriteParaview(def,tensionFields,fileName,1000,"_"); - - - // gsPiecewiseFunction<> membraneStresses; - // assembler->constructStress(mp_def,membraneStresses,stress_type::membrane); - // gsField<> membraneStress(mp_def,membraneStresses, true); - - // gsPiecewiseFunction<> flexuralStresses; - // assembler->constructStress(mp_def,flexuralStresses,stress_type::flexural); - // gsField<> flexuralStress(mp_def,flexuralStresses, true); - - // gsPiecewiseFunction<> stretches; - // assembler->constructStress(mp_def,stretches,stress_type::principal_stretch); - // gsField<> Stretches(mp_def,stretches, true); - - // gsPiecewiseFunction<> pstrain_m; - // assembler->constructStress(mp_def,pstrain_m,stress_type::principal_membrane_strain); - // gsField<> pstrainM(mp_def,pstrain_m, true); - - // gsPiecewiseFunction<> pstrain_f; - // assembler->constructStress(mp_def,pstrain_f,stress_type::principal_flexural_strain); - // gsField<> pstrainF(mp_def,pstrain_f, true); - - // gsPiecewiseFunction<> pstress_m; - // assembler->constructStress(mp_def,pstress_m,stress_type::principal_stress_membrane); - // gsField<> pstressM(mp_def,pstress_m, true); - - // gsPiecewiseFunction<> pstress_f; - // assembler->constructStress(mp_def,pstress_f,stress_type::principal_stress_flexural); - // gsField<> pstressF(mp_def,pstress_f, true); - - // gsPiecewiseFunction<> stretch1; - // assembler->constructStress(mp_def,stretch1,stress_type::principal_stretch_dir1); - // gsField<> stretchDir1(mp_def,stretch1, true); - - // gsPiecewiseFunction<> stretch2; - // assembler->constructStress(mp_def,stretch2,stress_type::principal_stretch_dir2); - // gsField<> stretchDir2(mp_def,stretch2, true); - - // gsPiecewiseFunction<> stretch3; - // assembler->constructStress(mp_def,stretch3,stress_type::principal_stretch_dir3); - // gsField<> stretchDir3(mp_def,stretch3, true); - - // gsPiecewiseFunction<> VMStresses; - // assembler->constructStress(mp_def,VMStresses,stress_type::von_mises_membrane); - // gsField<> VMStress(mp_def,VMStresses, true); - - - // gsWriteParaview(membraneStress,dirname + sep + "MembraneStress",5000); - // gsWriteParaview(VMStress,dirname + sep + "MembraneStressVM",5000); - // gsWriteParaview(Stretches,dirname + sep + "PrincipalStretch",5000); - // gsWriteParaview(pstrainM,dirname + sep + "PrincipalMembraneStrain",5000); - // gsWriteParaview(pstrainF,dirname + sep + "PrincipalFlexuralStrain",5000); - // gsWriteParaview(pstressM,dirname + sep + "PrincipalMembraneStress",5000); - // gsWriteParaview(pstressF,dirname + sep + "PrincipalFlexuralStress",5000); - // gsWriteParaview(stretchDir1,dirname + sep + "PrincipalDirection1",5000); - // gsWriteParaview(stretchDir1,dirname + sep + "PrincipalDirection1",5000); - // gsWriteParaview(stretchDir2,dirname + sep + "PrincipalDirection2",5000); - // gsWriteParaview(stretchDir3,dirname + sep + "PrincipalDirection3",5000); - - } - - delete assembler; - - return EXIT_SUCCESS; -} -#else//gsUnstructuredSplines_ENABLED -int main(int argc, char *argv[]) -{ - gsWarn<<"G+Smo is not compiled with the gsUnstructuredSplines module."; - return EXIT_FAILURE; -} -#endif -#else//gsKLShell_ENABLED -int main(int argc, char *argv[]) -{ - gsWarn<<"G+Smo is not compiled with the gsKLShell module."; - return EXIT_FAILURE; -} -#endif From 13789e50559418fddec69d3eb0f8da25460c39bd Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 25 Jun 2025 13:54:55 +0200 Subject: [PATCH 24/29] small review fix --- src/gsStaticSolvers/gsStaticDR.hpp | 2 +- src/gsStaticSolvers/gsStaticNewton.h | 4 ---- src/gsStaticSolvers/gsStaticNewton.hpp | 1 + 3 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/gsStaticSolvers/gsStaticDR.hpp b/src/gsStaticSolvers/gsStaticDR.hpp index f904038..26488b7 100644 --- a/src/gsStaticSolvers/gsStaticDR.hpp +++ b/src/gsStaticSolvers/gsStaticDR.hpp @@ -232,7 +232,7 @@ void gsStaticDR::_start() // m_residual = m_R.norm(); m_residual = m_forcing.norm(); // If the residual is 0 (e.g. with purely displacment loading), we set it to 1 to allow divisions - if (m_residual==0) m_residual=0; + if (m_residual==0) m_residual=1; // All residual norms are equal m_residualIni = m_residualOld = m_residual; } diff --git a/src/gsStaticSolvers/gsStaticNewton.h b/src/gsStaticSolvers/gsStaticNewton.h index e1fb04e..e99d78c 100644 --- a/src/gsStaticSolvers/gsStaticNewton.h +++ b/src/gsStaticSolvers/gsStaticNewton.h @@ -53,7 +53,6 @@ class gsStaticNewton : public gsStaticBase m_nonlinear(nullptr), m_residualFun(nullptr) { - m_dofs = m_linear.rows(); this->_init(); } @@ -77,7 +76,6 @@ class gsStaticNewton : public gsStaticBase m_residualFun(residual), m_ALresidualFun(nullptr) { - m_dofs = m_linear.rows(); m_dnonlinear = [this](gsVector const & x, gsVector const & /*dx*/, gsSparseMatrix & m) -> bool { return m_nonlinear(x,m); @@ -106,7 +104,6 @@ class gsStaticNewton : public gsStaticBase m_residualFun(nullptr), m_ALresidualFun(ALResidual) { - m_dofs = m_linear.rows(); m_L = 1.0; m_residualFun = [this](gsVector const & x, gsVector & result) -> bool { @@ -142,7 +139,6 @@ class gsStaticNewton : public gsStaticBase m_residualFun(residual), m_ALresidualFun(nullptr) { - m_dofs = m_linear.rows(); this->_init(); } diff --git a/src/gsStaticSolvers/gsStaticNewton.hpp b/src/gsStaticSolvers/gsStaticNewton.hpp index 0a47c6a..f03d4ca 100644 --- a/src/gsStaticSolvers/gsStaticNewton.hpp +++ b/src/gsStaticSolvers/gsStaticNewton.hpp @@ -259,6 +259,7 @@ void gsStaticNewton::_init() m_stabilityMethod = 0; m_start = false; + m_dofs = m_linear.rows(); if (m_dofs==0) gsWarn<<"The number of degrees of freedom is equal to zero. This can lead to bad initialization.\n"; From fa6926963009fa224b7b364ead2d186c189b9eea Mon Sep 17 00:00:00 2001 From: Hugo Verhelst <38376601+hverhelst@users.noreply.github.com> Date: Fri, 11 Jul 2025 09:53:58 +0200 Subject: [PATCH 25/29] Update benchmarks/benchmark_cylinder_DC.cpp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- benchmarks/benchmark_cylinder_DC.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 4492a42..6621c61 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -40,7 +40,7 @@ #include using namespace gismo; -void writeToCSVfile(std::string name, gsMatrix<> matrix) +void writeToCSVfile(const std::string& name, const gsMatrix<>& matrix) { std::ofstream file(name.c_str()); for(int i = 0; i < matrix.rows(); i++) From 30544937a33b9375831acae72e62317ab3e9b832 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst <38376601+hverhelst@users.noreply.github.com> Date: Fri, 11 Jul 2025 09:55:05 +0200 Subject: [PATCH 26/29] Update benchmarks/benchmark_cylinder_DC.cpp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- benchmarks/benchmark_cylinder_DC.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 6621c61..8f3a6ab 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -47,14 +47,13 @@ void writeToCSVfile(const std::string& name, const gsMatrix<>& matrix) { for(int j = 0; j < matrix.cols(); j++) { - std::string str = std::to_string(matrix(i,j)); if(j+1 == matrix.cols()) { - file< Date: Fri, 11 Jul 2025 10:14:46 +0200 Subject: [PATCH 27/29] add correct paths elasticity, fix frustrum example --- .../benchmark_Elasticity_Beam_APALM.cpp | 6 ++--- benchmarks/benchmark_FrustrumALM.cpp | 2 +- benchmarks/benchmark_Pillow.cpp | 2 +- benchmarks/benchmark_UniaxialTension.cpp | 26 ++++++++++++++++--- examples/gsElasticity_Modal.cpp | 6 ++--- examples/gsThinShell_Static_XML.cpp | 2 +- examples/snapping_element.cpp | 10 +++---- examples/snapping_element3D.cpp | 10 +++---- examples/snapping_example.cpp | 6 ++--- examples/snapping_example3D.cpp | 12 ++++----- filedata/pde/beam.xml | 8 ++++-- src/gsALMSolvers/gsALMCrisfield.hpp | 1 - tutorials/nonlinear_shell_arcLength.cpp | 2 +- tutorials/nonlinear_shell_dynamic.cpp | 2 +- tutorials/nonlinear_shell_static.cpp | 2 +- tutorials/nonlinear_solid_arcLength.cpp | 10 +++---- tutorials/nonlinear_solid_dynamic.cpp | 10 +++---- tutorials/nonlinear_solid_static.cpp | 10 +++---- 18 files changed, 75 insertions(+), 52 deletions(-) diff --git a/benchmarks/benchmark_Elasticity_Beam_APALM.cpp b/benchmarks/benchmark_Elasticity_Beam_APALM.cpp index d8dad99..1e3983d 100644 --- a/benchmarks/benchmark_Elasticity_Beam_APALM.cpp +++ b/benchmarks/benchmark_Elasticity_Beam_APALM.cpp @@ -14,9 +14,9 @@ #include #ifdef gsElasticity_ENABLED -#include -#include -#include +#include +#include +#include #endif #include diff --git a/benchmarks/benchmark_FrustrumALM.cpp b/benchmarks/benchmark_FrustrumALM.cpp index 7679e33..548aef8 100644 --- a/benchmarks/benchmark_FrustrumALM.cpp +++ b/benchmarks/benchmark_FrustrumALM.cpp @@ -110,7 +110,7 @@ int main (int argc, char** argv) // Initialise solution object mp_def = mp; - gsMultiBasis<> dbasis(mp,true); + gsMultiBasis<> dbasis(mp,false); gsInfo<<"Basis (patch 0): "<< mp.patch(0).basis() << "\n"; // Boundary conditions diff --git a/benchmarks/benchmark_Pillow.cpp b/benchmarks/benchmark_Pillow.cpp index f116d49..2983916 100644 --- a/benchmarks/benchmark_Pillow.cpp +++ b/benchmarks/benchmark_Pillow.cpp @@ -44,7 +44,7 @@ int main (int argc, char** argv) int verbose = 0; bool membrane = false; - real tolDR = 1e-1; + real_t tolDR = 1e-1; index_t Compressibility = 1; index_t material = 1; diff --git a/benchmarks/benchmark_UniaxialTension.cpp b/benchmarks/benchmark_UniaxialTension.cpp index c54dfd4..7e63078 100644 --- a/benchmarks/benchmark_UniaxialTension.cpp +++ b/benchmarks/benchmark_UniaxialTension.cpp @@ -3,7 +3,7 @@ @brief Uniaxial Tension Test benchmark e.g. Section 8.1. from Kiendl et al 2015 - Kiendl, J., Hsu, M.-C., Wu, M. C. H., & Reali, A. (2015). Isogeometric Kirchhoff–Love shell formulations for general hyperelastic materials. Computer Methods in Applied Mechanics and Engineering, 291, 280–303. https://doi.org/10.1016/J.CMA.2015.03.010 + Kiendl, J., Hsu, M.-C., Wu, M. C. H., & Reali, A. (2015). Isogeometric Kirchhoff–Love shell formlations for general hyperelastic materials. Computer Methods in Applied Mechanics and Engineering, 291, 280–303. https://doi.org/10.1016/J.CMA.2015.03.010 This file is part of the G+Smo library. @@ -344,11 +344,30 @@ int main (int argc, char** argv) assembler->constructSolution(solVector,mp_def); gsMatrix<> pts(2,1); - pts<<0.5,0.5; + pts<<1,0; lambdas.col(k) = assembler->computePrincipalStretches(pts,mp_def,0); S.at(k) = Lold / 1e-3 / lambdas(0,k) / lambdas(2,k); - San.at(k) = mu * (math::pow(lambdas(1,k),2)-1/lambdas(1,k)); + if (material==1 && !Compressibility) + San.at(k) = mu * (math::pow(lambdas(0,k),2)-1/lambdas(0,k)); + else if (material==3 && !Compressibility) + { + real_t c2 = 1.0 / (Ratio+1); + real_t c1 = 1.0 - c2; + San.at(k) =-mu*(c2*lambdas(0,k)*lambdas(0,k)+c2/lambdas(0,k)+c1)/lambdas(0,k)+lambdas(0,k)*(c1*lambdas(0,k)*mu+2*c2*mu); + } + else if (material==4 && !Compressibility) + { + real_t a1 = alpha1.value(0); + real_t a2 = alpha2.value(0); + real_t a3 = alpha3.value(0); + real_t m1 = mu1.value(0); + real_t m2 = mu2.value(0); + real_t m3 = mu3.value(0); + San.at(k) =-m1*math::pow((1./lambdas(0,k)),0.5*a1)-m2*math::pow((1./lambdas(0,k)),0.5*a2)-m3*math::pow((1./lambdas(0,k)),0.5*a3)+m1*math::pow(lambdas(0,k),a1)+m2*math::pow(lambdas(0,k),a2)+m3*math::pow(lambdas(0,k),a3); + } + else + San.at(k) = 0.; deformation = mp_def; deformation.patch(0).coefs() -= mp.patch(0).coefs();// assuming 1 patch here @@ -408,6 +427,7 @@ int main (int argc, char** argv) gsInfo<<"Lambdas:\n"< #ifdef gsElasticity_ENABLED -#include -#include -#include +#include +#include +#include #endif #include diff --git a/examples/gsThinShell_Static_XML.cpp b/examples/gsThinShell_Static_XML.cpp index 0576eb7..0c3b250 100644 --- a/examples/gsThinShell_Static_XML.cpp +++ b/examples/gsThinShell_Static_XML.cpp @@ -19,7 +19,7 @@ #endif #ifdef gsElasticity_ENABLED -#include +#include #endif #include diff --git a/examples/snapping_element.cpp b/examples/snapping_element.cpp index 9c2f490..12454ca 100644 --- a/examples/snapping_element.cpp +++ b/examples/snapping_element.cpp @@ -1,8 +1,8 @@ /** @file snapping_element.cpp @brief Arc-length analysis of a snapping meta material element. Inspired by - - Rafsanjani, A., Akbarzadeh, A., & Pasini, D. (2015). Snapping Mechanical Metamaterials under Tension. + + Rafsanjani, A., Akbarzadeh, A., & Pasini, D. (2015). Snapping Mechanical Metamaterials under Tension. Advanced Materials, 27(39), 5931–5935. https://doi.org/10.1002/adma.201502809 This file is part of the G+Smo library. @@ -22,9 +22,9 @@ #include #ifdef gsElasticity_ENABLED -#include -#include -#include +#include +#include +#include #endif using namespace gismo; diff --git a/examples/snapping_element3D.cpp b/examples/snapping_element3D.cpp index 89aaa5c..54d345e 100644 --- a/examples/snapping_element3D.cpp +++ b/examples/snapping_element3D.cpp @@ -1,8 +1,8 @@ /** @file snapping_element3D.cpp @brief Arc-length analysis of a snapping meta material element. Inspired by - - Rafsanjani, A., Akbarzadeh, A., & Pasini, D. (2015). Snapping Mechanical Metamaterials under Tension. + + Rafsanjani, A., Akbarzadeh, A., & Pasini, D. (2015). Snapping Mechanical Metamaterials under Tension. Advanced Materials, 27(39), 5931–5935. https://doi.org/10.1002/adma.201502809 This file is part of the G+Smo library. @@ -22,9 +22,9 @@ #include #ifdef gsElasticity_ENABLED -#include -#include -#include +#include +#include +#include #endif using namespace gismo; diff --git a/examples/snapping_example.cpp b/examples/snapping_example.cpp index 1e70c91..811d24d 100644 --- a/examples/snapping_example.cpp +++ b/examples/snapping_example.cpp @@ -22,9 +22,9 @@ #include #ifdef gsElasticity_ENABLED -#include -#include -#include +#include +#include +#include #endif using namespace gismo; diff --git a/examples/snapping_example3D.cpp b/examples/snapping_example3D.cpp index 375a185..ffd820a 100644 --- a/examples/snapping_example3D.cpp +++ b/examples/snapping_example3D.cpp @@ -1,8 +1,8 @@ /** @file snapping_example3D.cpp @brief Arc-length analysis of a snapping meta material. Inspired by - - Rafsanjani, A., Akbarzadeh, A., & Pasini, D. (2015). Snapping Mechanical Metamaterials under Tension. + + Rafsanjani, A., Akbarzadeh, A., & Pasini, D. (2015). Snapping Mechanical Metamaterials under Tension. Advanced Materials, 27(39), 5931–5935. https://doi.org/10.1002/adma.201502809 This file is part of the G+Smo library. @@ -22,9 +22,9 @@ #include #ifdef gsElasticity_ENABLED -#include -#include -#include +#include +#include +#include #endif using namespace gismo; @@ -332,7 +332,7 @@ int main(int argc, char *argv[]) }; assembler.options().setInt("MaterialLaw",material_law::neo_hooke_quad); - + //=============================================// // Solving // //=============================================// diff --git a/filedata/pde/beam.xml b/filedata/pde/beam.xml index f06216a..1c78e8f 100644 --- a/filedata/pde/beam.xml +++ b/filedata/pde/beam.xml @@ -38,7 +38,11 @@ - 0 + + 0 + 0 + 0 + 0 1 1 1 @@ -116,7 +120,7 @@ - + diff --git a/src/gsALMSolvers/gsALMCrisfield.hpp b/src/gsALMSolvers/gsALMCrisfield.hpp index 147d5e8..9af1454 100644 --- a/src/gsALMSolvers/gsALMCrisfield.hpp +++ b/src/gsALMSolvers/gsALMCrisfield.hpp @@ -341,7 +341,6 @@ void gsALMCrisfield::computeLambdas() m_note += "\tC"; // Compute eta computeLambdasModified(); - gsDebugVar(m_discriminant); if ((m_discriminant >= 0) && m_eta > 0.05) { // recompute lambdas with new eta diff --git a/tutorials/nonlinear_shell_arcLength.cpp b/tutorials/nonlinear_shell_arcLength.cpp index 8a0215d..9764eab 100644 --- a/tutorials/nonlinear_shell_arcLength.cpp +++ b/tutorials/nonlinear_shell_arcLength.cpp @@ -184,6 +184,6 @@ int main(int argc, char *argv[]) return EXIT_SUCCESS; #else GISMO_ERROR("The tutorial needs to be compiled with gsKLShell enabled"); - return EXIT_FAILED; + return EXIT_FAILURE; #endif }// end main \ No newline at end of file diff --git a/tutorials/nonlinear_shell_dynamic.cpp b/tutorials/nonlinear_shell_dynamic.cpp index d1af583..2602dd5 100644 --- a/tutorials/nonlinear_shell_dynamic.cpp +++ b/tutorials/nonlinear_shell_dynamic.cpp @@ -206,6 +206,6 @@ int main(int argc, char *argv[]) return EXIT_SUCCESS; #else GISMO_ERROR("The tutorial needs to be compiled with gsKLShell enabled"); - return EXIT_FAILED; + return EXIT_FAILURE; #endif }// end main \ No newline at end of file diff --git a/tutorials/nonlinear_shell_static.cpp b/tutorials/nonlinear_shell_static.cpp index ab6111c..1d2696c 100644 --- a/tutorials/nonlinear_shell_static.cpp +++ b/tutorials/nonlinear_shell_static.cpp @@ -194,6 +194,6 @@ int main(int argc, char *argv[]) return EXIT_SUCCESS; #else GISMO_ERROR("The tutorial needs to be compiled with gsKLShell enabled"); - return EXIT_FAILED; + return EXIT_FAILURE; #endif }// end main diff --git a/tutorials/nonlinear_solid_arcLength.cpp b/tutorials/nonlinear_solid_arcLength.cpp index 0119373..5fc6da9 100644 --- a/tutorials/nonlinear_solid_arcLength.cpp +++ b/tutorials/nonlinear_solid_arcLength.cpp @@ -15,9 +15,9 @@ //! [Includes] #ifdef gsElasticity_ENABLED -#include -#include -#include +#include +#include +#include #endif @@ -27,9 +27,9 @@ using namespace gismo; -#ifdef gsElasticity_ENABLED int main(int argc, char *argv[]) { +#ifdef gsElasticity_ENABLED //! [Parse command line] bool plot = false; index_t numRefine = 1; @@ -166,7 +166,7 @@ int main(int argc, char *argv[]) return EXIT_SUCCESS; #else GISMO_ERROR("The tutorial needs to be compiled with gsElasticity enabled"); - return EXIT_FAILED; + return EXIT_FAILURE; #endif }// end main \ No newline at end of file diff --git a/tutorials/nonlinear_solid_dynamic.cpp b/tutorials/nonlinear_solid_dynamic.cpp index cd9378a..3ba55b0 100644 --- a/tutorials/nonlinear_solid_dynamic.cpp +++ b/tutorials/nonlinear_solid_dynamic.cpp @@ -15,9 +15,9 @@ //! [Includes] #ifdef gsElasticity_ENABLED -#include -#include -#include +#include +#include +#include #endif @@ -27,9 +27,9 @@ using namespace gismo; -#ifdef gsElasticity_ENABLED int main(int argc, char *argv[]) { +#ifdef gsElasticity_ENABLED //! [Parse command line] bool plot = false; index_t numRefine = 1; @@ -193,7 +193,7 @@ for (index_t c=0; c!=3; c++) return EXIT_SUCCESS; #else GISMO_ERROR("The tutorial needs to be compiled with gsElasticity enabled"); - return EXIT_FAILED; + return EXIT_FAILURE; #endif }// end main \ No newline at end of file diff --git a/tutorials/nonlinear_solid_static.cpp b/tutorials/nonlinear_solid_static.cpp index 7529631..fa3f207 100644 --- a/tutorials/nonlinear_solid_static.cpp +++ b/tutorials/nonlinear_solid_static.cpp @@ -15,9 +15,9 @@ //! [Includes] #ifdef gsElasticity_ENABLED -#include -#include -#include +#include +#include +#include #endif @@ -27,9 +27,9 @@ using namespace gismo; -#ifdef gsElasticity_ENABLED int main(int argc, char *argv[]) { +#ifdef gsElasticity_ENABLED //! [Parse command line] bool plot = false; index_t numRefine = 1; @@ -175,7 +175,7 @@ int main(int argc, char *argv[]) return EXIT_SUCCESS; #else GISMO_ERROR("The tutorial needs to be compiled with gsElasticity enabled"); - return EXIT_FAILED; + return EXIT_FAILURE; #endif }// end main \ No newline at end of file From 9373594fba05e6f62202fa6bd474600ff1b91ac9 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Fri, 11 Jul 2025 11:08:25 +0200 Subject: [PATCH 28/29] small --- benchmarks/benchmark_cylinder_DC.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/benchmarks/benchmark_cylinder_DC.cpp b/benchmarks/benchmark_cylinder_DC.cpp index 8f3a6ab..699e2b8 100644 --- a/benchmarks/benchmark_cylinder_DC.cpp +++ b/benchmarks/benchmark_cylinder_DC.cpp @@ -1,7 +1,7 @@ -/** @file benchmark_TensionWrinkling.cpp +/** @file benchmark_cylinder_DC.cpp - @brief Computes the wrinkling behaviour of a thin sheet + @brief Computes the wrinkling behaviour of a cylinder or annulus This file is part of the G+Smo library. @@ -203,7 +203,7 @@ int main (int argc, char** argv) if (TFT) dirname = dirname + "_TFT"; if (membrane) - dirname = dirname + "_membrane"; + dirname = dirname + "_membrane"; gsFileManager::mkdir(dirname); From 7fa8dd91f7793c0de9860785c0bbcd710295daaa Mon Sep 17 00:00:00 2001 From: Hugo Verhelst <38376601+hverhelst@users.noreply.github.com> Date: Fri, 11 Jul 2025 11:09:02 +0200 Subject: [PATCH 29/29] Update benchmarks/benchmark_UniaxialTension.cpp --- benchmarks/benchmark_UniaxialTension.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/benchmark_UniaxialTension.cpp b/benchmarks/benchmark_UniaxialTension.cpp index 7e63078..98d16dc 100644 --- a/benchmarks/benchmark_UniaxialTension.cpp +++ b/benchmarks/benchmark_UniaxialTension.cpp @@ -3,7 +3,7 @@ @brief Uniaxial Tension Test benchmark e.g. Section 8.1. from Kiendl et al 2015 - Kiendl, J., Hsu, M.-C., Wu, M. C. H., & Reali, A. (2015). Isogeometric Kirchhoff–Love shell formlations for general hyperelastic materials. Computer Methods in Applied Mechanics and Engineering, 291, 280–303. https://doi.org/10.1016/J.CMA.2015.03.010 +Kiendl, J., Hsu, M.-C., Wu, M. C. H., & Reali, A. (2015). Isogeometric Kirchhoff–Love shell formulations for general hyperelastic materials. Computer Methods in Applied Mechanics and Engineering, 291, 280–303. https://doi.org/10.1016/J.CMA.2015.03.010 This file is part of the G+Smo library.