Diffusion akantu (features/30-implementing-a-checkpoint-restart-mechanism) Browse

30-implementin…

structural_mechanics_model.cc
View Options

			/**
			* @file structural_mechanics_model.cc
			*
			* @author Nicolas Richart <nicolas.richart@epfl.ch>
			* @author Damien Spielmann <damien.spielmann@epfl.ch>
			* @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
			* @author Fabian Barras <fabian.barras@epfl.ch>
			*
			* @date creation: Fri Jul 15 2011
			* @date last modification: Mon Sep 15 2014
			*
			* @brief Model implementation for StucturalMechanics elements
			*
			* @section LICENSE
			*
			* Copyright (©) 2014 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 <http://www.gnu.org/licenses/>.
			*
			*/

			/* -------------------------------------------------------------------------- */
			#include "structural_mechanics_model.hh"
			#include "group_manager_inline_impl.cc"
			#include "dumpable_inline_impl.hh"
			#include "aka_math.hh"
			#include "integration_scheme_2nd_order.hh"
			#include "static_communicator.hh"
			#include "sparse_matrix.hh"
			#include "solver.hh"

			#ifdef AKANTU_USE_MUMPS
			#include "solver_mumps.hh"
			#endif

			#ifdef AKANTU_USE_IOHELPER
			# include "dumper_paraview.hh"
			# include "dumper_elemental_field.hh"
			#endif
			/* -------------------------------------------------------------------------- */

			__BEGIN_AKANTU__

			const StructuralMechanicsModelOptions default_structural_mechanics_model_options(_static);

			/* -------------------------------------------------------------------------- */
			StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh,
			UInt dim,
			const ID & id,
			const MemoryID & memory_id) :
			Model(mesh, dim, id, memory_id),
			time_step(NAN), f_m2a(1.0),
			stress("stress", id, memory_id),
			element_material("element_material", id, memory_id),
			set_ID("beam sets", id, memory_id),
			stiffness_matrix(NULL),
			mass_matrix(NULL),
			velocity_damping_matrix(NULL),
			jacobian_matrix(NULL),
			solver(NULL),
			rotation_matrix("rotation_matices", id, memory_id),
			increment_flag(false),
			integrator(NULL) {
			AKANTU_DEBUG_IN();

			registerFEEngineObject<MyFEEngineType>("StructuralMechanicsFEEngine", mesh, spatial_dimension);

			this->displacement_rotation = NULL;
			this->mass = NULL;
			this->velocity = NULL;
			this->acceleration = NULL;
			this->force_momentum = NULL;
			this->residual = NULL;
			this->blocked_dofs = NULL;
			this->increment = NULL;

			this->previous_displacement = NULL;

			if(spatial_dimension == 2)
			nb_degree_of_freedom = 3;
			else if (spatial_dimension == 3)
			nb_degree_of_freedom = 6;
			else {
			AKANTU_DEBUG_TO_IMPLEMENT();
			}

			this->mesh.registerDumper<DumperParaview>("paraview_all", id, true);
			this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			StructuralMechanicsModel::~StructuralMechanicsModel() {
			AKANTU_DEBUG_IN();

			delete integrator;
			delete solver;
			delete stiffness_matrix;
			delete jacobian_matrix;
			delete mass_matrix;
			delete velocity_damping_matrix;

			AKANTU_DEBUG_OUT();
			}

			void StructuralMechanicsModel::setTimeStep(Real time_step) {
			this->time_step = time_step;

			#if defined(AKANTU_USE_IOHELPER)
			this->mesh.getDumper().setTimeStep(time_step);
			#endif
			}

			/* -------------------------------------------------------------------------- */
			/* Initialisation */
			/* -------------------------------------------------------------------------- */

			void StructuralMechanicsModel::initFull(const ModelOptions & options) {

			const StructuralMechanicsModelOptions & stmm_options =
			dynamic_cast<const StructuralMechanicsModelOptions &>(options);

			method = stmm_options.analysis_method;

			initModel();
			initArrays();

			displacement_rotation->clear();
			velocity ->clear();
			acceleration ->clear();
			force_momentum ->clear();
			residual ->clear();
			blocked_dofs ->clear();
			increment ->clear();


			Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural);
			Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural);
			for (; it != end; ++it) {
			computeRotationMatrix(*it);
			}

			switch(method) {
			case _implicit_dynamic:
			initImplicit();
			break;
			case _static:
			initSolver();
			break;
			default:
			AKANTU_EXCEPTION("analysis method not recognised by StructuralMechanicsModel");
			break;
			}
			}

			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::initArrays() {
			AKANTU_DEBUG_IN();

			UInt nb_nodes = getFEEngine().getMesh().getNbNodes();
			std::stringstream sstr_disp; sstr_disp << id << ":displacement";
			std::stringstream sstr_velo; sstr_velo << id << ":velocity";
			std::stringstream sstr_acce; sstr_acce << id << ":acceleration";
			std::stringstream sstr_forc; sstr_forc << id << ":force";
			std::stringstream sstr_resi; sstr_resi << id << ":residual";
			std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs";
			std::stringstream sstr_incr; sstr_incr << id << ":increment";

			displacement_rotation = &(alloc<Real>(sstr_disp.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
			velocity = &(alloc<Real>(sstr_velo.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
			acceleration = &(alloc<Real>(sstr_acce.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
			force_momentum = &(alloc<Real>(sstr_forc.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
			residual = &(alloc<Real>(sstr_resi.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));
			blocked_dofs = &(alloc<bool>(sstr_boun.str(), nb_nodes, nb_degree_of_freedom, false));
			increment = &(alloc<Real>(sstr_incr.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE));

			Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural);
			Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural);
			for (; it != end; ++it) {
			UInt nb_element = getFEEngine().getMesh().getNbElement(*it);
			UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(*it);

			element_material.alloc(nb_element, 1, *it, _not_ghost);
			set_ID.alloc(nb_element, 1, *it, _not_ghost);

			UInt size = getTangentStiffnessVoigtSize(*it);
			stress.alloc(nb_element * nb_quadrature_points, size , *it, _not_ghost);
			}

			dof_synchronizer = new DOFSynchronizer(getFEEngine().getMesh(), nb_degree_of_freedom);
			dof_synchronizer->initLocalDOFEquationNumbers();
			dof_synchronizer->initGlobalDOFEquationNumbers();

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::initModel() {
			getFEEngine().initShapeFunctions(_not_ghost);
			}

			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) {
			AKANTU_DEBUG_IN();

			const Mesh & mesh = getFEEngine().getMesh();
			#if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake
			AKANTU_DEBUG_ERROR("You should at least activate one solver.");
			#else
			UInt nb_global_node = mesh.getNbGlobalNodes();

			std::stringstream sstr; sstr << id << ":jacobian_matrix";
			jacobian_matrix = new SparseMatrix(nb_global_node * nb_degree_of_freedom, _symmetric, sstr.str(), memory_id);

			jacobian_matrix->buildProfile(mesh, *dof_synchronizer, nb_degree_of_freedom);

			std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix";
			stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id);

			#ifdef AKANTU_USE_MUMPS
			std::stringstream sstr_solv; sstr_solv << id << ":solver";
			solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());

			dof_synchronizer->initScatterGatherCommunicationScheme();
			#else
			AKANTU_DEBUG_ERROR("You should at least activate one solver.");
			#endif //AKANTU_USE_MUMPS

			solver->initialize(options);
			#endif //AKANTU_HAS_SOLVER

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) {
			AKANTU_DEBUG_IN();

			if (!increment) setIncrementFlagOn();

			initSolver(solver_options);

			if(integrator) delete integrator;
			integrator = new TrapezoidalRule2();

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			UInt StructuralMechanicsModel::getTangentStiffnessVoigtSize(const ElementType & type) {
			UInt size;
			#define GET_TANGENT_STIFFNESS_VOIGT_SIZE(type) \
			size = getTangentStiffnessVoigtSize<type>();

			AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_TANGENT_STIFFNESS_VOIGT_SIZE);
			#undef GET_TANGENT_STIFFNESS_VOIGT_SIZE

			return size;
			}

			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::assembleStiffnessMatrix() {
			AKANTU_DEBUG_IN();

			stiffness_matrix->clear();

			Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural);
			Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural);
			for (; it != end; ++it) {
			ElementType type = *it;

			#define ASSEMBLE_STIFFNESS_MATRIX(type) \
			assembleStiffnessMatrix<type>();

			AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX);
			#undef ASSEMBLE_STIFFNESS_MATRIX
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template<>
			void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_2>(Array<Real> & rotations){

			ElementType type = _bernoulli_beam_2;
			Mesh & mesh = getFEEngine().getMesh();
			UInt nb_element = mesh.getNbElement(type);

			Array<UInt>::iterator< Vector<UInt> > connec_it = mesh.getConnectivity(type).begin(2);
			Array<Real>::vector_iterator nodes_it = mesh.getNodes().begin(spatial_dimension);
			Array<Real>::matrix_iterator R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom);

			for (UInt e = 0; e < nb_element; ++e, ++R_it, ++connec_it) {
			Matrix<Real> & R = *R_it;
			Vector<UInt> & connec = *connec_it;

			Vector<Real> x2;
			x2 = nodes_it[connec(1)]; // X2
			Vector<Real> x1;
			x1 = nodes_it[connec(0)]; // X1

			Real le = x1.distance(x2);
			Real c = (x2(0) - x1(0)) / le;
			Real s = (x2(1) - x1(1)) / le;

			/// Definition of the rotation matrix
			R(0,0) = c; R(0,1) = s; R(0,2) = 0.;
			R(1,0) = -s; R(1,1) = c; R(1,2) = 0.;
			R(2,0) = 0.; R(2,1) = 0.; R(2,2) = 1.;
			}
			}


			/* -------------------------------------------------------------------------- */
			template<>
			void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_3>(Array<Real> & rotations){
			ElementType type = _bernoulli_beam_3;
			Mesh & mesh = getFEEngine().getMesh();
			UInt nb_element = mesh.getNbElement(type);

			Array<Real>::vector_iterator n_it = mesh.getNormals(type).begin(spatial_dimension);
			Array<UInt>::iterator< Vector<UInt> > connec_it = mesh.getConnectivity(type).begin(2);
			Array<Real>::vector_iterator nodes_it = mesh.getNodes().begin(spatial_dimension);

			Matrix<Real> Pe (spatial_dimension, spatial_dimension);
			Matrix<Real> Pg (spatial_dimension, spatial_dimension);
			Matrix<Real> inv_Pg(spatial_dimension, spatial_dimension);
			Vector<Real> x_n(spatial_dimension); // x vect n

			Array<Real>::matrix_iterator R_it =
			rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom);

			for (UInt e=0 ; e < nb_element; ++e, ++n_it, ++connec_it, ++R_it) {
			Vector<Real> & n = *n_it;
			Matrix<Real> & R = *R_it;
			Vector<UInt> & connec = *connec_it;

			Vector<Real> x;
			x = nodes_it[connec(1)]; // X2
			Vector<Real> y;
			y = nodes_it[connec(0)]; // X1

			Real l = x.distance(y);
			x -= y; // X2 - X1
			x_n.crossProduct(x, n);

			Pe.eye();
			Pe(0, 0) *= l;
			Pe(1, 1) *= -l;

			Pg(0,0) = x(0); Pg(0,1) = x_n(0); Pg(0,2) = n(0);
			Pg(1,0) = x(1); Pg(1,1) = x_n(1); Pg(1,2) = n(1);
			Pg(2,0) = x(2); Pg(2,1) = x_n(2); Pg(2,2) = n(2);

			inv_Pg.inverse(Pg);

			Pe *= inv_Pg;
			for (UInt i = 0; i < spatial_dimension; ++i) {
			for (UInt j = 0; j < spatial_dimension; ++j) {
			R(i, j) = Pe(i, j);
			R(i + spatial_dimension,j + spatial_dimension) = Pe(i, j);
			}
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			template<>
			void StructuralMechanicsModel::computeRotationMatrix<_kirchhoff_shell>(Array<Real> & rotations){
			ElementType type = _kirchhoff_shell;
			Mesh & mesh = getFEEngine().getMesh();
			UInt nb_element = mesh.getNbElement(type);

			Array<UInt>::iterator< Vector<UInt> > connec_it = mesh.getConnectivity(type).begin(3);
			Array<Real>::vector_iterator nodes_it = mesh.getNodes().begin(spatial_dimension);

			Matrix<Real> Pe (spatial_dimension, spatial_dimension);
			Matrix<Real> Pg (spatial_dimension, spatial_dimension);
			Matrix<Real> inv_Pg(spatial_dimension, spatial_dimension);

			Array<Real>::matrix_iterator R_it =
			rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom);

			for (UInt e=0 ; e < nb_element; ++e, ++connec_it, ++R_it) {

			Pe.eye();

			Matrix<Real> & R = *R_it;
			Vector<UInt> & connec = *connec_it;

			Vector<Real> x2;
			x2 = nodes_it[connec(1)]; // X2
			Vector<Real> x1;
			x1 = nodes_it[connec(0)]; // X1
			Vector<Real> x3;
			x3 = nodes_it[connec(2)]; // X3




			Vector<Real> Pg_col_1=x2-x1;

			Vector<Real> Pg_col_2 = x3-x1;

			Vector<Real> Pg_col_3(spatial_dimension);
			Pg_col_3.crossProduct(Pg_col_1, Pg_col_2);


			for (UInt i = 0; i <spatial_dimension; ++i) {
			Pg(i,0) = Pg_col_1(i);
			Pg(i,1) = Pg_col_2(i);
			Pg(i,2) = Pg_col_3(i);
			}


			inv_Pg.inverse(Pg);
			// Pe *= inv_Pg;
			Pe.eye();

			for (UInt i = 0; i < spatial_dimension; ++i) {
			for (UInt j = 0; j < spatial_dimension; ++j) {
			R(i, j) = Pe(i, j);
			R(i + spatial_dimension,j + spatial_dimension) = Pe(i, j);
			}
			}

			}
			}

			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::computeRotationMatrix(const ElementType & type) {
			Mesh & mesh = getFEEngine().getMesh();

			UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
			UInt nb_element = mesh.getNbElement(type);

			if(!rotation_matrix.exists(type)) {
			rotation_matrix.alloc(nb_element,
			nb_degree_of_freedomnb_nodes_per_element nb_degree_of_freedom*nb_nodes_per_element,
			type);
			} else {
			rotation_matrix(type).resize(nb_element);
			}
			rotation_matrix(type).clear();

			Array<Real>rotations(nb_element, nb_degree_of_freedom * nb_degree_of_freedom);
			rotations.clear();

			#define COMPUTE_ROTATION_MATRIX(type) \
			computeRotationMatrix<type>(rotations);

			AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX);
			#undef COMPUTE_ROTATION_MATRIX


			Array<Real>::matrix_iterator R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom);
			Array<Real>::matrix_iterator T_it =
			rotation_matrix(type).begin(nb_degree_of_freedom*nb_nodes_per_element,
			nb_degree_of_freedom*nb_nodes_per_element);

			for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) {
			Matrix<Real> & T = *T_it;
			Matrix<Real> & R = *R_it;
			T.clear();
			for (UInt k = 0; k < nb_nodes_per_element; ++k){
			for (UInt i = 0; i < nb_degree_of_freedom; ++i)
			for (UInt j = 0; j < nb_degree_of_freedom; ++j)
			T(knb_degree_of_freedom + i, knb_degree_of_freedom + j) = R(i, j);
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::computeStresses() {
			AKANTU_DEBUG_IN();

			Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural);
			Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural);
			for (; it != end; ++it) {
			ElementType type = *it;

			#define COMPUTE_STRESS_ON_QUAD(type) \
			computeStressOnQuad<type>();

			AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD);
			#undef COMPUTE_STRESS_ON_QUAD
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::updateResidual() {
			AKANTU_DEBUG_IN();
			residual->copy(*force_momentum);

			Array<Real> ku(*displacement_rotation, true);
			ku = stiffness_matrix;
			*residual -= ku;

			this->updateResidualInternal();

			AKANTU_DEBUG_OUT();

			}

			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::implicitPred() {
			AKANTU_DEBUG_IN();

			if(previous_displacement) previous_displacement->copy(*displacement_rotation);

			if(method == _implicit_dynamic)
			integrator->integrationSchemePred(time_step,
			*displacement_rotation,
			*velocity,
			*acceleration,
			*blocked_dofs);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::implicitCorr() {
			AKANTU_DEBUG_IN();

			Real * incr_val = increment->storage();
			bool * boun_val = blocked_dofs->storage();

			UInt nb_nodes = displacement_rotation->getSize();

			for (UInt j = 0; j < nb_nodes * nb_degree_of_freedom; ++j, ++incr_val, ++boun_val){
			incr_val = (1. - *boun_val);
			}

			if(method == _implicit_dynamic) {
			integrator->integrationSchemeCorrDispl(time_step,
			*displacement_rotation,
			*velocity,
			*acceleration,
			*blocked_dofs,
			*increment);
			} else {

			Real * disp_val = displacement_rotation->storage();
			incr_val = increment->storage();

			for (UInt j = 0; j < nb_nodes *nb_degree_of_freedom; ++j, ++disp_val, ++incr_val){
			disp_val += incr_val;
			}
			}

			AKANTU_DEBUG_OUT();
			}
			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::updateResidualInternal() {
			AKANTU_DEBUG_IN();

			AKANTU_DEBUG_INFO("Update the residual");
			// f = f_ext - f_int - Ma - Cv = r - Ma - Cv;

			if(method != _static) {
			// f -= Ma
			if(mass_matrix) {
			// if full mass_matrix
			Array<Real> * Ma = new Array<Real>(*acceleration, true, "Ma");
			Ma = *mass_matrix;
			/// \todo check unit conversion for implicit dynamics
			// *Ma /= f_m2a
			residual -= Ma;
			delete Ma;
			} else if (mass) {

			// else lumped mass
			UInt nb_nodes = acceleration->getSize();
			UInt nb_degree_of_freedom = acceleration->getNbComponent();

			Real * mass_val = mass->storage();
			Real * accel_val = acceleration->storage();
			Real * res_val = residual->storage();
			bool * blocked_dofs_val = blocked_dofs->storage();

			for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
			if(!(*blocked_dofs_val)) {
			res_val -= accel_val * *mass_val /f_m2a;
			}
			blocked_dofs_val++;
			res_val++;
			mass_val++;
			accel_val++;
			}
			} else {
			AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
			}

			// f -= Cv
			if(velocity_damping_matrix) {
			Array<Real> * Cv = new Array<Real>(*velocity);
			Cv = *velocity_damping_matrix;
			residual -= Cv;
			delete Cv;
			}
			}

			AKANTU_DEBUG_OUT();
			}


			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::setIncrementFlagOn() {
			AKANTU_DEBUG_IN();

			if(!increment) {
			UInt nb_nodes = mesh.getNbNodes();
			std::stringstream sstr_inc; sstr_inc << id << ":increment";
			increment = &(alloc<Real>(sstr_inc.str(), nb_nodes, spatial_dimension, 0.));
			}

			increment_flag = true;

			AKANTU_DEBUG_OUT();
			}


			/* -------------------------------------------------------------------------- */
			void StructuralMechanicsModel::solve() {
			AKANTU_DEBUG_IN();

			AKANTU_DEBUG_INFO("Solving an implicit step.");

			UInt nb_nodes = displacement_rotation->getSize();

			/// todo residual = force - Kxr * d_bloq
			jacobian_matrix->copyContent(*stiffness_matrix);
			jacobian_matrix->applyBoundary(*blocked_dofs);

			jacobian_matrix->saveMatrix("Kbound.mtx");

			increment->clear();

			solver->setRHS(*residual);

			solver->factorize();
			solver->solve(*increment);

			Real * increment_val = increment->storage();
			Real * displacement_val = displacement_rotation->storage();
			bool * blocked_dofs_val = blocked_dofs->storage();

			for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
			if(!(*blocked_dofs_val)) {
			displacement_val += increment_val;
			}
			else {
			*increment_val = 0.0;
			}

			displacement_val++;
			blocked_dofs_val++;
			increment_val++;
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template<>
			bool StructuralMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){
			AKANTU_DEBUG_IN();

			UInt nb_nodes = displacement_rotation->getSize();
			UInt nb_degree_of_freedom = displacement_rotation->getNbComponent();

			error = 0;
			Real norm[2] = {0., 0.};
			Real * increment_val = increment->storage();
			bool * blocked_dofs_val = blocked_dofs->storage();
			Real * displacement_val = displacement_rotation->storage();

			for (UInt n = 0; n < nb_nodes; ++n) {
			for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
			if(!(*blocked_dofs_val)) {
			norm[0] += increment_val *increment_val;
			norm[1] += displacement_val *displacement_val;
			}
			blocked_dofs_val++;
			increment_val++;
			displacement_val++;
			}
			}

			StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum);

			norm[0] = sqrt(norm[0]);
			norm[1] = sqrt(norm[1]);
			AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something went wrong in the solve phase");

			AKANTU_DEBUG_OUT();
			if(norm[1] > Math::getTolerance())
			error = norm[0] / norm[1];
			else
			error = norm[0]; //In case the total displacement is zero!
			return (error < tolerance);
			}

			/* -------------------------------------------------------------------------- */
			template<>
			bool StructuralMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) {
			AKANTU_DEBUG_IN();

			UInt nb_nodes = residual->getSize();

			norm = 0;
			Real * residual_val = residual->storage();
			bool * blocked_dofs_val = blocked_dofs->storage();

			for (UInt n = 0; n < nb_nodes; ++n) {
			bool is_local_node = mesh.isLocalOrMasterNode(n);
			if(is_local_node) {
			for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
			if(!(*blocked_dofs_val)) {
			norm += residual_val *residual_val;
			}
			blocked_dofs_val++;
			residual_val++;
			}
			} else {
			blocked_dofs_val += spatial_dimension;
			residual_val += spatial_dimension;
			}
			}

			StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum);

			norm = sqrt(norm);

			AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase");

			AKANTU_DEBUG_OUT();
			return (norm < tolerance);
			}


			/* -------------------------------------------------------------------------- */
			bool StructuralMechanicsModel::testConvergenceIncrement(Real tolerance) {
			Real error;
			bool tmp = testConvergenceIncrement(tolerance, error);

			AKANTU_DEBUG_INFO("Norm of increment : " << error);

			return tmp;
			}

			/* -------------------------------------------------------------------------- */
			bool StructuralMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error) {
			AKANTU_DEBUG_IN();
			Mesh & mesh= getFEEngine().getMesh();
			UInt nb_nodes = displacement_rotation->getSize();
			UInt nb_degree_of_freedom = displacement_rotation->getNbComponent();

			Real norm = 0;
			Real * increment_val = increment->storage();
			bool * blocked_dofs_val = blocked_dofs->storage();

			for (UInt n = 0; n < nb_nodes; ++n) {
			bool is_local_node = mesh.isLocalOrMasterNode(n);
			for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
			if(!(*blocked_dofs_val) && is_local_node) {
			norm += increment_val *increment_val;
			}
			blocked_dofs_val++;
			increment_val++;
			}
			}

			StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum);

			error = sqrt(norm);
			AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase");

			AKANTU_DEBUG_OUT();
			return (error < tolerance);
			}

			/* -------------------------------------------------------------------------- */
			template<>
			void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_2>(Array<Real> & tangent_moduli) {
			UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_2);
			UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_2);
			UInt tangent_size = 2;

			Array<Real>::matrix_iterator D_it = tangent_moduli.begin(tangent_size, tangent_size);
			Array<UInt> & el_mat = element_material(_bernoulli_beam_2, _not_ghost);

			for (UInt e = 0; e < nb_element; ++e) {
			UInt mat = el_mat(e);
			Real E = materials[mat].E;
			Real A = materials[mat].A;
			Real I = materials[mat].I;
			for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) {
			Matrix<Real> & D = *D_it;
			D(0,0) = E * A;
			D(1,1) = E * I;
			}
			}
			}
			/* -------------------------------------------------------------------------- */
			template<>
			void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_3>(Array<Real> & tangent_moduli) {
			UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_3);
			UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_3);
			UInt tangent_size = 4;

			Array<Real>::matrix_iterator D_it = tangent_moduli.begin(tangent_size, tangent_size);

			for (UInt e = 0; e < nb_element; ++e) {
			UInt mat = element_material(_bernoulli_beam_3, _not_ghost)(e);
			Real E = materials[mat].E;
			Real A = materials[mat].A;
			Real Iz = materials[mat].Iz;
			Real Iy = materials[mat].Iy;
			Real GJ = materials[mat].GJ;
			for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) {
			Matrix<Real> & D = *D_it;
			D(0,0) = E * A;
			D(1,1) = E * Iz;
			D(2,2) = E * Iy;
			D(3,3) = GJ;
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			template<>
			void StructuralMechanicsModel::computeTangentModuli<_kirchhoff_shell>(Array<Real> & tangent_moduli) {
			UInt nb_element = getFEEngine().getMesh().getNbElement(_kirchhoff_shell);
			UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_kirchhoff_shell);
			UInt tangent_size = 6;

			Array<Real>::matrix_iterator D_it = tangent_moduli.begin(tangent_size, tangent_size);

			for (UInt e = 0; e < nb_element; ++e) {
			UInt mat = element_material(_kirchhoff_shell, _not_ghost)(e);
			Real E = materials[mat].E;
			Real nu = materials[mat].nu;
			Real t = materials[mat].t;

			Real HH=Et/(1-nunu);
			Real DD=Ettt/(12(1-nu*nu));

			for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) {
			Matrix<Real> & D = *D_it;
			D(0,0) = HH;
			D(0,1) = HH*nu;
			D(1,0) = HH*nu;
			D(1,1) = HH;
			D(2,2) = HH*(1-nu)/2;
			D(3,3) = DD;
			D(3,4) = DD*nu;
			D(4,3) = DD*nu;
			D(4,4) = DD;
			D(5,5) = DD*(1-nu)/2;
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			template<>
			void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_bernoulli_beam_2>(Array<Real> & b, bool local) {
			UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_2);
			UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_bernoulli_beam_2);
			UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_2);

			MyFEEngineType & fem = getFEEngineClass<MyFEEngineType>();
			Array<Real>::const_vector_iterator shape_Np = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 0).begin(nb_nodes_per_element);
			Array<Real>::const_vector_iterator shape_Mpp = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 1).begin(nb_nodes_per_element);
			Array<Real>::const_vector_iterator shape_Lpp = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 2).begin(nb_nodes_per_element);

			UInt tangent_size = getTangentStiffnessVoigtSize<_bernoulli_beam_2>();
			UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom;
			b.clear();
			Array<Real>::matrix_iterator B_it = b.begin(tangent_size, bt_d_b_size);

			for (UInt e = 0; e < nb_element; ++e) {
			for (UInt q = 0; q < nb_quadrature_points; ++q) {
			Matrix<Real> & B = *B_it;
			const Vector<Real> & Np = *shape_Np;
			const Vector<Real> & Lpp = *shape_Lpp;
			const Vector<Real> & Mpp = *shape_Mpp;

			B(0,0) = Np(0);
			B(0,3) = Np(1);

			B(1,1) = Mpp(0);
			B(1,2) = Lpp(0);
			B(1,4) = Mpp(1);
			B(1,5) = Lpp(1);

			++B_it;
			++shape_Np;
			++shape_Mpp;
			++shape_Lpp;
			}

			// ++R_it;
			}
			}

			/* -------------------------------------------------------------------------- */
			template<>
			void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_bernoulli_beam_3>(Array<Real> & b, bool local) {
			MyFEEngineType & fem = getFEEngineClass<MyFEEngineType>();

			UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_3);
			UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_bernoulli_beam_3);
			UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_3);

			Array<Real>::const_vector_iterator shape_Np = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 0).begin(nb_nodes_per_element);
			Array<Real>::const_vector_iterator shape_Mpp = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 1).begin(nb_nodes_per_element);
			Array<Real>::const_vector_iterator shape_Lpp = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 2).begin(nb_nodes_per_element);

			UInt tangent_size = getTangentStiffnessVoigtSize<_bernoulli_beam_3>();
			UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom;

			b.clear();

			Array<Real>::matrix_iterator B_it = b.begin(tangent_size, bt_d_b_size);

			for (UInt e = 0; e < nb_element; ++e) {
			for (UInt q = 0; q < nb_quadrature_points; ++q) {
			Matrix<Real> & B = *B_it;

			const Vector<Real> & Np = *shape_Np;
			const Vector<Real> & Lpp = *shape_Lpp;
			const Vector<Real> & Mpp = *shape_Mpp;

			B(0,0) = Np(0);
			B(0,6) = Np(1);

			B(1,1) = Mpp(0);
			B(1,5) = Lpp(0);
			B(1,7) = Mpp(1);
			B(1,11) = Lpp(1);

			B(2,2) = Mpp(0);
			B(2,4) = -Lpp(0);
			B(2,8) = Mpp(1);
			B(2,10) = -Lpp(1);

			B(3,3) = Np(0);
			B(3,9) = Np(1);

			++B_it;
			++shape_Np;
			++shape_Mpp;
			++shape_Lpp;
			}
			}
			}
			/* -------------------------------------------------------------------------- */
			template<>
			void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_kirchhoff_shell>(Array<Real> & b, bool local) {
			MyFEEngineType & fem = getFEEngineClass<MyFEEngineType>();


			UInt nb_element = getFEEngine().getMesh().getNbElement(_kirchhoff_shell);
			UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_kirchhoff_shell);
			UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_kirchhoff_shell);



			Array<Real>::const_matrix_iterator shape_Np = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 0).begin(2,nb_nodes_per_element);
			Array<Real>::const_matrix_iterator shape_Nx1p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 1).begin(2,nb_nodes_per_element);
			Array<Real>::const_matrix_iterator shape_Nx2p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 2).begin(2,nb_nodes_per_element);
			Array<Real>::const_matrix_iterator shape_Nx3p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 3).begin(2,nb_nodes_per_element);
			Array<Real>::const_matrix_iterator shape_Ny1p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 4).begin(2,nb_nodes_per_element);
			Array<Real>::const_matrix_iterator shape_Ny2p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 5).begin(2,nb_nodes_per_element);
			Array<Real>::const_matrix_iterator shape_Ny3p = fem.getShapesDerivatives(_kirchhoff_shell, _not_ghost, 6).begin(2,nb_nodes_per_element);

			UInt tangent_size = getTangentStiffnessVoigtSize<_kirchhoff_shell>();
			UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom;

			b.clear();

			Array<Real>::matrix_iterator B_it = b.begin(tangent_size, bt_d_b_size);

			for (UInt e = 0; e < nb_element; ++e) {
			for (UInt q = 0; q < nb_quadrature_points; ++q) {
			Matrix<Real> & B = *B_it;

			const Matrix<Real> & Np = *shape_Np;
			const Matrix<Real> & Nx1p = *shape_Nx1p;
			const Matrix<Real> & Nx2p = *shape_Nx2p;
			const Matrix<Real> & Nx3p = *shape_Nx3p;
			const Matrix<Real> & Ny1p = *shape_Ny1p;
			const Matrix<Real> & Ny2p = *shape_Ny2p;
			const Matrix<Real> & Ny3p = *shape_Ny3p;

			B(0,0) = Np(0,0);
			B(0,6) = Np(0,1);
			B(0,12) = Np(0,2);

			B(1,1) = Np(1,0);
			B(1,7) = Np(1,1);
			B(1,13) = Np(1,2);

			B(2,0) = Np(1,0);
			B(2,1) = Np(0,0);
			B(2,6) = Np(1,1);
			B(2,7) = Np(0,1);
			B(2,12) = Np(1,2);
			B(2,13) = Np(0,2);

			B(3,2) = Nx1p(0,0);
			B(3,3) = Nx2p(0,0);
			B(3,4) = Nx3p(0,0);
			B(3,8) = Nx1p(0,1);
			B(3,9) = Nx2p(0,1);
			B(3,10) = Nx3p(0,1);
			B(3,14) = Nx1p(0,2);
			B(3,15) = Nx2p(0,2);
			B(3,16) = Nx3p(0,2);

			B(4,2) = Ny1p(1,0);
			B(4,3) = Ny2p(1,0);
			B(4,4) = Ny3p(1,0);
			B(4,8) = Ny1p(1,1);
			B(4,9) = Ny2p(1,1);
			B(4,10) = Ny3p(1,1);
			B(4,14) = Ny1p(1,2);
			B(4,15) = Ny2p(1,2);
			B(4,16) = Ny3p(1,2);

			B(5,2) = Nx1p(1,0) + Ny1p(0,0);
			B(5,3) = Nx2p(1,0) + Ny2p(0,0);
			B(5,4) = Nx3p(1,0) + Ny3p(0,0);
			B(5,8) = Nx1p(1,1) + Ny1p(0,1);
			B(5,9) = Nx2p(1,1) + Ny2p(0,1);
			B(5,10) = Nx3p(1,1) + Ny3p(0,1);
			B(5,14) = Nx1p(1,2) + Ny1p(0,2);
			B(5,15) = Nx2p(1,2) + Ny2p(0,2);
			B(5,16) = Nx3p(1,2) + Ny3p(0,2);

			++B_it;
			++shape_Np;
			++shape_Nx1p;
			++shape_Nx2p;
			++shape_Nx3p;
			++shape_Ny1p;
			++shape_Ny2p;
			++shape_Ny3p;
			}
			}
			}


			/* -------------------------------------------------------------------------- */

			dumper::Field * StructuralMechanicsModel
			::createNodalFieldBool(const std::string & field_name,
			const std::string & group_name,
			bool padding_flag) {

			std::map<std::string,Array<bool>* > uint_nodal_fields;
			uint_nodal_fields["blocked_dofs" ] = blocked_dofs;

			dumper::Field * field = NULL;
			field = mesh.createNodalField(uint_nodal_fields[field_name],group_name);
			return field;

			}

			/* -------------------------------------------------------------------------- */


			dumper::Field * StructuralMechanicsModel
			::createNodalFieldReal(const std::string & field_name,
			const std::string & group_name,
			bool padding_flag) {

			UInt n;
			if(spatial_dimension == 2) {
			n = 2;
			} else n = 3;

			dumper::Field * field = NULL;

			if (field_name == "displacement"){
			field = mesh.createStridedNodalField(displacement_rotation,group_name,n,0,padding_flag);
			}
			if (field_name == "rotation"){
			field = mesh.createStridedNodalField(displacement_rotation,group_name,
			nb_degree_of_freedom-n,n,padding_flag);
			}
			if (field_name == "force"){
			field = mesh.createStridedNodalField(force_momentum,group_name,n,0,padding_flag);
			}
			if (field_name == "momentum"){
			field = mesh.createStridedNodalField(force_momentum,group_name,
			nb_degree_of_freedom-n,n,padding_flag);
			}
			if (field_name == "residual"){
			field = mesh.createNodalField(residual,group_name,padding_flag);
			}

			return field;

			}

			/* -------------------------------------------------------------------------- */

			dumper::Field * StructuralMechanicsModel
			::createElementalField(const std::string & field_name,
			const std::string & group_name,
			bool padding_flag,
			const ElementKind & kind,
			const std::string & fe_engine_id){


			dumper::Field * field = NULL;

			if(field_name == "element_index_by_material")
			field = mesh.createElementalField<UInt, Vector, dumper::ElementalField >(field_name,group_name,this->spatial_dimension,kind);

			return field;
			}

			/* -------------------------------------------------------------------------- */





			__END_AKANTU__

rAKA/src/model/structural_mechanics/structural_mechanics_model.cc9054f18b641efeatures/30-implementin…

structural_mechanics_model.ccView Options

rAKA/src/model/structural_mechanics/structural_mechanics_model.cc
9054f18b641efeatures/30-implementin…

structural_mechanics_model.cc
View Options