diff --git a/src/model/structural_mechanics/structural_mechanics_model.hh b/src/model/structural_mechanics/structural_mechanics_model.hh index 600705cb9..659223e19 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.hh +++ b/src/model/structural_mechanics/structural_mechanics_model.hh @@ -1,391 +1,392 @@ /** * @file structural_mechanics_model.hh * * @author Fabian Barras * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Wed Apr 22 2015 * * @brief Particular implementation of the structural elements in the *StructuralMechanicsModel * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__ #define __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "model.hh" #include "integrator_gauss.hh" #include "shape_linked.hh" #include "aka_types.hh" #include "dumpable.hh" #include "solver.hh" #include "integration_scheme_2nd_order.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class SparseMatrix; } __BEGIN_AKANTU__ struct StructuralMaterial { Real E; Real A; Real I; Real Iz; Real Iy; Real GJ; Real rho; Real t; Real nu; }; struct StructuralMechanicsModelOptions : public ModelOptions { StructuralMechanicsModelOptions(AnalysisMethod analysis_method = _static) : analysis_method(analysis_method) {} AnalysisMethod analysis_method; }; extern const StructuralMechanicsModelOptions default_structural_mechanics_model_options; class StructuralMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEEngineTemplate MyFEEngineType; StructuralMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "structural_mechanics_model", const MemoryID & memory_id = 0); virtual ~StructuralMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize fully the model void initFull(const ModelOptions & options = default_structural_mechanics_model_options); /// initialize the internal vectors void initArrays(); /// initialize the model void initModel(); /// initialize the solver void initSolver(SolverOptions & options = _solver_no_options); /// initialize the stuff for the implicit solver void initImplicit(bool dynamic = false, SolverOptions & solver_options = _solver_no_options); /// compute the stresses per elements void computeStresses(); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); /// implicit time integration predictor void implicitPred(); /// implicit time integration corrector void implicitCorr(); /// update the residual vector void updateResidual(); /// solve the system void solve(); bool testConvergenceIncrement(Real tolerance); bool testConvergenceIncrement(Real tolerance, Real & error); void computeRotationMatrix(const ElementType & type); protected: UInt getTangentStiffnessVoigtSize(const ElementType & type); /// compute Rotation Matrices template void computeRotationMatrix(__attribute__((unused)) Array & rotations){}; /// compute A and solve @f[ A\delta u = f_ext - f_int @f] template - void solve(Array & increment, Real block_val = 1.); + void solve(Array & increment, + __attribute__((unused)) Real block_val = 1.); /* ------------------------------------------------------------------------ */ /* Mass (structural_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// computes rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// finish the computation of residual to solve in increment void updateResidualInternal(); /* ------------------------------------------------------------------------ */ private: template inline UInt getTangentStiffnessVoigtSize(); template void assembleStiffnessMatrix(); template void assembleMass(); template void computeStressOnQuad(); template void computeTangentModuli(Array & tangent_moduli); template void transferBMatrixToSymVoigtBMatrix(Array & B, bool local = false); template void transferNMatrixToSymVoigtNMatrix(Array & N_matrix); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: virtual dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag); virtual dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag); virtual dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const ElementKind & kind, const std::string & fe_engine_id = ""); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step void setTimeStep(Real time_step); /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the StructuralMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement_rotation, Array &); /// get the StructuralMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the StructuralMechanicsModel::acceleration vector, updated /// by /// StructuralMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the StructuralMechanicsModel::force vector (boundary forces) AKANTU_GET_MACRO(Force, *force_momentum, Array &); /// get the StructuralMechanicsModel::residual vector, computed by /// StructuralMechanicsModel::updateResidual AKANTU_GET_MACRO(Residual, *residual, const Array &); /// get the StructuralMechanicsModel::boundary vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(RotationMatrix, rotation_matrix, Real); AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, const SparseMatrix &); AKANTU_GET_MACRO(MassMatrix, *mass_matrix, const SparseMatrix &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Set_ID, set_ID, UInt); void addMaterial(StructuralMaterial & material) { materials.push_back(material); } /** * @brief set the StructuralMechanicsModel::increment_flag to on, the * activate the * update of the StructuralMechanicsModel::increment vector */ void setIncrementFlagOn(); /* ------------------------------------------------------------------------ */ /* Boundaries (structural_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: /// Compute Linear load function set in global axis template void computeForcesByGlobalTractionArray(const Array & tractions); /// Compute Linear load function set in local axis template void computeForcesByLocalTractionArray(const Array & tractions); /// compute force vector from a function(x,y,momentum) that describe stresses template void computeForcesFromFunction(BoundaryFunction in_function, BoundaryFunctionType function_type); /** * solve a step (predictor + convergence loop + corrector) using the * the given convergence method (see akantu::SolveConvergenceMethod) * and the given convergence criteria (see * akantu::SolveConvergenceCriteria) **/ template bool solveStep(Real tolerance, UInt max_iteration = 100); template bool solveStep(Real tolerance, Real & error, UInt max_iteration = 100); /// test if the system is converged template bool testConvergence(Real tolerance, Real & error); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement_rotation; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement; /// velocities array Array * velocity; /// accelerations array Array * acceleration; /// forces array Array * force_momentum; /// lumped mass array Array * mass; /// stress arraz ElementTypeMapArray stress; /// residuals array Array * residual; /// boundaries array Array * blocked_dofs; /// position of a dof in the K matrix Array * equation_number; ElementTypeMapArray element_material; // Define sets of beams ElementTypeMapArray set_ID; /// local equation_number to global unordered_map::type local_eq_num_to_global; /// stiffness matrix SparseMatrix * stiffness_matrix; /// mass matrix SparseMatrix * mass_matrix; /// velocity damping matrix SparseMatrix * velocity_damping_matrix; /// jacobian matrix SparseMatrix * jacobian_matrix; /// increment of displacement Array * increment; /// solver for implicit Solver * solver; /// number of degre of freedom UInt nb_degree_of_freedom; // Rotation matrix ElementTypeMapArray rotation_matrix; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// flag defining if the increment must be computed or not bool increment_flag; /// integration scheme of second order used IntegrationScheme2ndOrder * integrator; /* -------------------------------------------------------------------------- */ std::vector materials; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model_inline_impl.cc" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const StructuralMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__ */ diff --git a/src/solver/static_solver.cc b/src/solver/static_solver.cc index a76686b79..d606dc902 100644 --- a/src/solver/static_solver.cc +++ b/src/solver/static_solver.cc @@ -1,124 +1,125 @@ /** * @file static_solver.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Mon Oct 13 2014 * @date last modification: Fri Oct 16 2015 * * @brief implementation of the static solver * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "static_solver.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_PETSC #include #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ StaticSolver::StaticSolver() : CommunicatorEventHandler(), is_initialized(false) { StaticCommunicator::getStaticCommunicator().registerEventHandler(*this); } /* -------------------------------------------------------------------------- */ StaticSolver::~StaticSolver() { --this->nb_references; if (this->nb_references == 0) { StaticCommunicator::getStaticCommunicator().unregisterEventHandler(*this); delete this->static_solver; } } /* -------------------------------------------------------------------------- */ StaticSolver & StaticSolver::getStaticSolver() { if (nb_references == 0) static_solver = new StaticSolver(); ++nb_references; return *static_solver; } #ifdef AKANTU_USE_PETSC #if PETSC_VERSION_MAJOR >= 3 && PETSC_VERSION_MINOR >= 5 static PetscErrorCode PETScErrorHandler(MPI_Comm, int line, const char * dir, const char * file, PetscErrorCode number, PetscErrorType type, const char * message, void *) { AKANTU_DEBUG_ERROR("An error occured in PETSc in file \"" << file << ":" << line << "\" - PetscErrorCode " << number << " - \"" << message << "\""); } #else static PetscErrorCode PETScErrorHandler(MPI_Comm, int line, __attribute__((unused)) const char * func, __attribute__((unused)) const char * dir, const char * file, PetscErrorCode number, __attribute__((unused)) PetscErrorType type, const char * message, void *) { AKANTU_DEBUG_ERROR("An error occured in PETSc in file \"" << file << ":" << line << "\" - PetscErrorCode " << number << " - \"" << message << "\""); } #endif #endif /* -------------------------------------------------------------------------- */ -void StaticSolver::initialize(int & argc, char **& argv) { +void StaticSolver::initialize(__attribute__((unused)) int & argc, + __attribute__((unused)) char **& argv) { if (this->is_initialized) return; // AKANTU_DEBUG_ASSERT(this->is_initialized != true, "The static solver has // already been initialized"); #ifdef AKANTU_USE_PETSC PetscErrorCode petsc_error = PetscInitialize(&argc, &argv, NULL, NULL); if (petsc_error != 0) { AKANTU_DEBUG_ERROR( "An error occured while initializing Petsc (PetscErrorCode " << petsc_error << ")"); } PetscPushErrorHandler(PETScErrorHandler, NULL); #endif this->is_initialized = true; } /* -------------------------------------------------------------------------- */ void StaticSolver::finalize() { ParentEventHandler::sendEvent( StaticSolverEvent::BeforeStaticSolverDestroyEvent()); AKANTU_DEBUG_ASSERT(this->is_initialized == true, "The static solver has not been initialized"); #ifdef AKANTU_USE_PETSC PetscFinalize(); #endif this->is_initialized = false; } __END_AKANTU__ diff --git a/test/test_model/test_non_local_toolbox/test_material.hh b/test/test_model/test_non_local_toolbox/test_material.hh index 04d3fa9b9..e6883d764 100644 --- a/test/test_model/test_non_local_toolbox/test_material.hh +++ b/test/test_model/test_non_local_toolbox/test_material.hh @@ -1,78 +1,73 @@ /** * @file test_material.hh * * @author Aurelia Isabel Cuba Ramos * * @date creation: Thu Feb 21 2013 * @date last modification: Wed Nov 25 2015 * * @brief test material for the non-local neighborhood base test * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_damage.hh" #include "material_non_local.hh" #ifndef __TEST_MATERIAL_HH__ #define __TEST_MATERIAL_HH__ __BEGIN_AKANTU__ -template +template class TestMaterial : public MaterialDamage, - public MaterialNonLocal{ - -/* -------------------------------------------------------------------------- */ -/* Constructor/Destructor */ -/* -------------------------------------------------------------------------- */ - + public MaterialNonLocal { + /* ------------------------------------------------------------------------ */ + /* Constructor/Destructor */ + /* ------------------------------------------------------------------------ */ public: - TestMaterial(SolidMechanicsModel & model, const ID & id); - virtual ~TestMaterial() {}; + virtual ~TestMaterial(){}; typedef MaterialNonLocal MyNonLocalParent; typedef MaterialDamage MyLocalParent; -/* -------------------------------------------------------------------------- */ -/* Methods */ -/* -------------------------------------------------------------------------- */ + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ public: void initMaterial(); - void computeNonLocalStresses(GhostType ghost_type) {}; + void computeNonLocalStresses(__attribute__((unused)) GhostType ghost_type){}; void insertQuadsInNeighborhoods(GhostType ghost_type); virtual void registerNeighborhood(); -protected: - -/* -------------------------------------------------------------------------- */ -/* Members */ -/* -------------------------------------------------------------------------- */ + /* ------------------------------------------------------------------------ */ + /* Members */ + /* ------------------------------------------------------------------------ */ private: - InternalField grad_u_nl; - + InternalField grad_u_nl; }; __END_AKANTU__ #endif /* __TEST_MATERIAL_HH__ */ diff --git a/test/test_model/test_non_local_toolbox/test_material_damage.hh b/test/test_model/test_non_local_toolbox/test_material_damage.hh index 2cec37f8d..01d51bdf3 100644 --- a/test/test_model/test_non_local_toolbox/test_material_damage.hh +++ b/test/test_model/test_non_local_toolbox/test_material_damage.hh @@ -1,79 +1,72 @@ /** * @file test_material_damage.hh * * @author Aurelia Isabel Cuba Ramos * * @date creation: Thu Feb 21 2013 * @date last modification: Thu Oct 15 2015 * * @brief test material damage for the non-local remove damage test * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_damage.hh" #include "material_non_local.hh" #ifndef __TEST_MATERIAL_DAMAGE_HH__ #define __TEST_MATERIAL_DAMAGE_HH__ __BEGIN_AKANTU__ -template -class TestMaterialDamage : public MaterialDamage, - public MaterialNonLocal { - -/* -------------------------------------------------------------------------- */ -/* Constructor/Destructor */ -/* -------------------------------------------------------------------------- */ - +template +class TestMaterialDamage : public MaterialDamage, + public MaterialNonLocal { + /* ------------------------------------------------------------------------ */ + /* Constructor/Destructor */ + /* ------------------------------------------------------------------------ */ public: - TestMaterialDamage(SolidMechanicsModel & model, const ID & id); - virtual ~TestMaterialDamage() {}; + virtual ~TestMaterialDamage(){}; typedef MaterialNonLocal MyNonLocalParent; -/* -------------------------------------------------------------------------- */ -/* Methods */ -/* -------------------------------------------------------------------------- */ + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ public: void initMaterial(); - //void computeNonLocalStress(ElementType type, GhostType ghost_type = _not_ghost); - - void computeNonLocalStresses(GhostType ghost_type) {}; + // void computeNonLocalStress(ElementType type, GhostType ghost_type = + // _not_ghost); + void computeNonLocalStresses(__attribute__((unused)) GhostType ghost_type){}; void insertQuadsInNeighborhoods(GhostType ghost_type); -protected: - - -/* -------------------------------------------------------------------------- */ -/* Members */ -/* -------------------------------------------------------------------------- */ + /* ------------------------------------------------------------------------ */ + /* Members */ + /* ------------------------------------------------------------------------ */ private: - InternalField grad_u_nl; - + InternalField grad_u_nl; }; __END_AKANTU__ #endif /* __TEST_MATERIAL_DAMAGE_HH__ */ diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc index 484e6e148..1c60226cc 100644 --- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc @@ -1,195 +1,206 @@ /** * @file test_structural_mechanics_model_bernoulli_beam_dynamics.cc * * @author Sébastien Hartmann * * @date creation: Mon Jul 07 2014 * @date last modification: Sun Oct 19 2014 * * @brief Test for _bernouilli_beam in dynamic * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh_struct.hh" #include "structural_mechanics_model.hh" #include "material.hh" #include using namespace akantu; /* -------------------------------------------------------------------------- */ #define TYPE _bernoulli_beam_2 - -static Real analytical_solution(Real time, Real L, Real rho, Real E, Real A, Real I, Real F) { - Real omega = M_PI*M_PI/L/L*sqrt(E*I/rho); +static Real analytical_solution(Real time, Real L, Real rho, Real E, + __attribute__((unused)) Real A, Real I, + Real F) { + Real omega = M_PI * M_PI / L / L * sqrt(E * I / rho); Real sum = 0.; - UInt i=5; - for(UInt n=1; n <= i; n +=2) { - sum += (1. - cos(n*n*omega*time))/pow(n, 4); + UInt i = 5; + for (UInt n = 1; n <= i; n += 2) { + sum += (1. - cos(n * n * omega * time)) / pow(n, 4); } - return 2.*F*pow(L, 3)/pow(M_PI, 4)/E/I*sum; + return 2. * F * pow(L, 3) / pow(M_PI, 4) / E / I * sum; } - -//load +// load const Real F = 0.5e4; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ -int main(int argc, char *argv[]){ +int main(int argc, char * argv[]) { initialize(argc, argv); Mesh beams(2); debug::setDebugLevel(dblWarning); const ElementType type = _bernoulli_beam_2; - /* -------------------------------------------------------------------------- */ + /* -------------------------------------------------------------------------- + */ // Mesh UInt nb_element = 8; UInt nb_nodes = nb_element + 1; Real total_length = 10.; - Real length = total_length/nb_element; + Real length = total_length / nb_element; Real heigth = 1.; - + Array & nodes = const_cast &>(beams.getNodes()); nodes.resize(nb_nodes); - + beams.addConnectivityType(type); - Array &connectivity = const_cast &>(beams.getConnectivity(type)); + Array & connectivity = + const_cast &>(beams.getConnectivity(type)); connectivity.resize(nb_element); beams.initNormals(); Array & normals = const_cast &>(beams.getNormals(type)); normals.resize(nb_element); - - for (UInt i = 0; i < nb_nodes; ++i){ - nodes(i,0) = i * length; - nodes(i,1) = 0; + + for (UInt i = 0; i < nb_nodes; ++i) { + nodes(i, 0) = i * length; + nodes(i, 1) = 0; } - - for (UInt i = 0; i < nb_element; ++i){ - connectivity(i,0) = i; - connectivity(i,1) = i+1; + + for (UInt i = 0; i < nb_element; ++i) { + connectivity(i, 0) = i; + connectivity(i, 1) = i + 1; } - - /* -------------------------------------------------------------------------- */ + + /* -------------------------------------------------------------------------- + */ // Materials - + StructuralMechanicsModel model(beams); - + StructuralMaterial mat1; mat1.E = 120e6; mat1.rho = 1000; mat1.A = heigth; - mat1.I = heigth*heigth*heigth/12; + mat1.I = heigth * heigth * heigth / 12; model.addMaterial(mat1); - /* -------------------------------------------------------------------------- */ + /* -------------------------------------------------------------------------- + */ // Forces // model.initFull(); model.initFull(StructuralMechanicsModelOptions(_implicit_dynamic)); - const Array &position = beams.getNodes(); - Array &forces = model.getForce(); - Array &displacement = model.getDisplacement(); - Array &boundary = model.getBlockedDOFs(); + const Array & position = beams.getNodes(); + Array & forces = model.getForce(); + Array & displacement = model.getDisplacement(); + Array & boundary = model.getBlockedDOFs(); UInt node_to_print = -1; forces.clear(); displacement.clear(); // boundary.clear(); - //model.getElementMaterial(type)(i,0) = 0; - //model.getElementMaterial(type)(i,0) = 1; + // model.getElementMaterial(type)(i,0) = 0; + // model.getElementMaterial(type)(i,0) = 1; for (UInt i = 0; i < nb_element; ++i) { - model.getElementMaterial(type)(i,0) = 0; + model.getElementMaterial(type)(i, 0) = 0; } - for (UInt n =0; n(load, _bft_traction) - /* -------------------------------------------------------------------------- */ + /* -------------------------------------------------------------------------- + */ // Boundary conditions // true ~ displacement is blocked - boundary(0,0) = true; - boundary(0,1) = true; - boundary(nb_nodes-1,1) = true; - /* -------------------------------------------------------------------------- */ + boundary(0, 0) = true; + boundary(0, 1) = true; + boundary(nb_nodes - 1, 1) = true; + /* -------------------------------------------------------------------------- + */ // "Solve" - + Real time = 0; model.assembleStiffnessMatrix(); model.assembleMass(); model.dump(); model.getStiffnessMatrix().saveMatrix("K.mtx"); model.getMassMatrix().saveMatrix("M.mt"); Real time_step = 1e-4; model.setTimeStep(time_step); - std::cout << "Time" << " | " << "Mid-Span Displacement" << std::endl; + std::cout << "Time" + << " | " + << "Mid-Span Displacement" << std::endl; /// time loop for (UInt s = 1; time < 0.64; ++s) { model.solveStep<_scm_newton_raphson_tangent, _scc_increment>(1e-12, 1000); - - pos << s << "," << time << "," << displacement(node_to_print, 1) << "," << analytical_solution(s*time_step, total_length, mat1.rho, mat1.E, mat1.A, mat1.I, F) << std::endl; - // pos << s << "," << time << "," << displacement(node_to_print, 1) << "," << analytical_solution(s*time_step) << std::endl; + + pos << s << "," << time << "," << displacement(node_to_print, 1) << "," + << analytical_solution(s * time_step, total_length, mat1.rho, mat1.E, + mat1.A, mat1.I, F) << std::endl; + // pos << s << "," << time << "," << displacement(node_to_print, 1) << + // "," << analytical_solution(s*time_step) << std::endl; time += time_step; - if(s % 100 == 0) - std::cout << time << " | " << displacement(node_to_print, 1) << std::endl; + if (s % 100 == 0) + std::cout << time << " | " << displacement(node_to_print, 1) + << std::endl; model.dump(); } pos.close(); finalize(); return EXIT_SUCCESS; } -