diff --git a/cmake/Modules/FindMumps.cmake b/cmake/Modules/FindMumps.cmake index 14e4f3d98..0db04fb96 100644 --- a/cmake/Modules/FindMumps.cmake +++ b/cmake/Modules/FindMumps.cmake @@ -1,293 +1,315 @@ #=============================================================================== # @file FindMumps.cmake # # @author Nicolas Richart # # @date creation: Fri Oct 24 2014 # @date last modification: Wed Jan 13 2016 # # @brief The find_package file for the Mumps 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 . # #=============================================================================== set(_MUMPS_COMPONENTS "sequential" "parallel" "double" "float" "complex_double" "complex_float") if(NOT Mumps_FIND_COMPONENTS) set(Mumps_FIND_COMPONENTS "parallel" "double" "float" "complex_double" "complex_float") endif() #=============================================================================== enable_language(Fortran) set(MUMPS_PRECISIONS) set(MUMPS_PLAT) foreach(_comp ${Mumps_FIND_COMPONENTS}) if("${_comp}" STREQUAL "sequential") set(MUMPS_PLAT _seq) #default plat on debian based distribution endif() if("${_comp}" STREQUAL "float") list(APPEND MUMPS_PRECISIONS s) endif() if("${_comp}" STREQUAL "double") list(APPEND MUMPS_PRECISIONS d) endif() if("${_comp}" STREQUAL "complex_float") list(APPEND MUMPS_PRECISIONS c) endif() if("${_comp}" STREQUAL "complex_double") list(APPEND MUMPS_PRECISIONS z) endif() endforeach() if(NOT MUMPS_PRECISIONS) set(MUMPS_PRECISIONS s d c z) endif() list(GET MUMPS_PRECISIONS 0 _first_precision) string(TOUPPER "${_first_precision}" _u_first_precision) find_path(MUMPS_INCLUDE_DIR ${_first_precision}mumps_c.h PATHS "${MUMPS_DIR}" ENV MUMPS_DIR PATH_SUFFIXES include ) mark_as_advanced(MUMPS_INCLUDE_DIR) set(_mumps_required_vars) foreach(_precision ${MUMPS_PRECISIONS}) string(TOUPPER "${_precision}" _u_precision) find_library(MUMPS_LIBRARY_${_u_precision}MUMPS NAMES ${_precision}mumps${MUMPS_PREFIX} PATHS "${MUMPS_DIR}" ENV MUMPS_DIR PATH_SUFFIXES lib ) mark_as_advanced(MUMPS_LIBRARY_${_u_precision}MUMPS) list(APPEND _mumps_required_vars MUMPS_LIBRARY_${_u_precision}MUMPS) list(APPEND MUMPS_LIBRARIES_ALL ${MUMPS_LIBRARY_${_u_precision}MUMPS}) endforeach() if(MUMPS_LIBRARY_${_u_first_precision}MUMPS MATCHES ".*${_first_precision}mumps.*${CMAKE_STATIC_LIBRARY_SUFFIX}") # Assuming mumps was compiled as a static library set(MUMPS_LIBRARY_TYPE STATIC CACHE INTERNAL "" FORCE) if (CMAKE_Fortran_COMPILER MATCHES ".*gfortran") set(_compiler_specific gfortran) elseif (CMAKE_Fortran_COMPILER MATCHES ".*ifort") set(_compiler_specific ifcore) else() message("Compiler ${CMAKE_Fortran_COMPILER} is not known, you will probably " "have to add semething instead of this message to be able to test mumps " "install") endif() else() set(MUMPS_LIBRARY_TYPE SHARED CACHE INTERNAL "" FORCE) endif() function(mumps_add_dependency _pdep _libs) string(TOUPPER ${_pdep} _u_pdep) if(_pdep STREQUAL "mumps_common") find_library(MUMPS_LIBRARY_COMMON mumps_common${MUMPS_PREFIX} PATHS "${MUMPS_DIR}" ENV MUMPS_DIR PATH_SUFFIXES lib ) set(${_libs} ${MUMPS_LIBRARY_COMMON} PARENT_SCOPE) mark_as_advanced(MUMPS_LIBRARY_COMMON) elseif(_pdep STREQUAL "pord") find_library(MUMPS_LIBRARY_PORD pord${MUMPS_PREFIX} PATHS "${MUMPS_DIR}" ENV MUMPS_DIR PATH_SUFFIXES lib ) set(${_libs} ${MUMPS_LIBRARY_PORD} PARENT_SCOPE) mark_as_advanced(MUMPS_LIBRARY_PORD) elseif(_pdep MATCHES "Scotch") find_package(Scotch REQUIRED ${ARGN} QUIET) if(ARGN) list(GET ARGN 1 _comp) string(TOUPPER ${_comp} _u_comp) set(${_libs} ${SCOTCH_LIBRARY_${_u_comp}} PARENT_SCOPE) else() set(${_libs} ${${_u_pdep}_LIBRARIES} PARENT_SCOPE) endif() elseif(_pdep MATCHES "MPI") if(MUMPS_PLAT STREQUAL "_seq") find_library(MUMPS_LIBRARY_MPISEQ mpiseq${MUMPS_PREFIX} PATHS "${MUMPS_DIR}" ENV MUMPS_DIR PATH_SUFFIXES lib ) + if (NOT MUMPS_LIBRARY_MPISEQ) + message("Could not find libmpiseq for sequential version of MUMPS, was " + "MUMPS compiled in sequential ?") + endif() set(${_libs} ${MUMPS_LIBRARY_MPISEQ} PARENT_SCOPE) mark_as_advanced(MUMPS_LIBRARY_MPISEQ) else() find_package(MPI REQUIRED C Fortran QUIET) set(${_libs} ${MPI_C_LIBRARIES} ${MPI_Fortran_LIBRARIES} PARENT_SCOPE) endif() elseif(_pdep MATCHES "Threads") find_package(Threads REQUIRED) set(${_libs} Threads::Threads PARENT_SCOPE) + elseif(_pdep MATCHES "OpenMP") + find_package(OpenMP REQUIRED) + set(${_libs} OpenMP::OpenMP_C PARENT_SCOPE) + elseif(_pdep MATCHES "Math") + set(${_libs} m PARENT_SCOPE) else() find_package(${_pdep} REQUIRED QUIET) set(${_libs} ${${_u_pdep}_LIBRARIES} ${${_u_pdep}_LIBRARY} PARENT_SCOPE) endif() endfunction() function(mumps_find_dependencies) - set(_libraries_all ${MUMPS_LIBRARIES_ALL}) + set(_libraries_all m ${MUMPS_LIBRARIES_ALL}) set(_include_dirs ${MUMPS_INCLUDE_DIR}) set(_mumps_test_dir "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}") file(READ ${CMAKE_CURRENT_LIST_DIR}/CheckFindMumps.c _output) file(WRITE "${_mumps_test_dir}/mumps_test_code.c" "#include <${_first_precision}mumps_c.h> ${_u_first_precision}MUMPS_STRUC_C id; #define mumps_c ${_first_precision}mumps_c #define Real ${_u_first_precision}MUMPS_REAL ") if(MUMPS_PLAT STREQUAL _seq) file(APPEND "${_mumps_test_dir}/mumps_test_code.c" "#define MUMPS_SEQ ") else() file(APPEND "${_mumps_test_dir}/mumps_test_code.c" "// #undef MUMPS_SEQ ") find_package(MPI REQUIRED) list(APPEND _compiler_specific ${MPI_C_LIBRARIES}) list(APPEND _include_dirs ${MPI_C_INCLUDE_PATH} ${MPI_INCLUDE_DIR}) endif() file(APPEND "${_mumps_test_dir}/mumps_test_code.c" "${_output}") #=============================================================================== set(_mumps_dep_symbol_BLAS ${_first_precision}gemm) set(_mumps_dep_symbol_ScaLAPACK numroc) + set(_mumps_dep_symbol_LAPACK ilaenv) set(_mumps_dep_symbol_MPI mpi_send) set(_mumps_dep_symbol_Scotch SCOTCH_graphInit) set(_mumps_dep_symbol_Scotch_ptscotch scotchfdgraphexit) set(_mumps_dep_symbol_Scotch_esmumps esmumps) set(_mumps_dep_symbol_mumps_common mumps_abort) set(_mumps_dep_symbol_pord SPACE_ordering) set(_mumps_dep_symbol_METIS metis_nodend) set(_mumps_dep_symbol_Threads pthread_create) + set(_mumps_dep_symbol_OpenMP GOMP_loop_end_nowait) + # TODO find missing symbols for IOMP + set(_mumps_dep_symbol_Math lround) set(_mumps_dep_symbol_ParMETIS ParMETIS_V3_NodeND) # added for fucking macosx that cannot fail at link set(_mumps_run_dep_symbol_mumps_common mumps_fac_descband) set(_mumps_run_dep_symbol_MPI mpi_bcast) set(_mumps_run_dep_symbol_ScaLAPACK idamax) set(_mumps_dep_comp_Scotch_ptscotch COMPONENTS ptscotch) set(_mumps_dep_comp_Scotch_esmumps COMPONENTS esmumps) - set(_mumps_potential_dependencies mumps_common pord BLAS ScaLAPACK MPI - Scotch Scotch_ptscotch Scotch_esmumps METIS ParMETIS Threads) + set(_mumps_potential_dependencies + mumps_common pord + BLAS LAPACK ScaLAPACK + MPI + Scotch Scotch_ptscotch Scotch_esmumps + METIS ParMETIS + Threads OpenMP + Math) #=============================================================================== set(_retry_try_run TRUE) set(_retry_count 0) # trying only as long as we add dependencies to avoid inifinte loop in case of an unkown dependency while (_retry_try_run AND _retry_count LESS 100) - try_run(_mumps_run _mumps_compiles "${_mumps_test_dir}" "${_mumps_test_dir}/mumps_test_code.c" + try_run(_mumps_run _mumps_compiles + "${_mumps_test_dir}" + "${_mumps_test_dir}/mumps_test_code.c" CMAKE_FLAGS "-DINCLUDE_DIRECTORIES:STRING=${_include_dirs}" LINK_LIBRARIES ${_libraries_all} ${_libraries_all} ${_compiler_specific} RUN_OUTPUT_VARIABLE _run COMPILE_OUTPUT_VARIABLE _out) set(_retry_compile FALSE) #message("COMPILATION outputs: \n${_out} \n RUN OUTPUT \n${_run}") if(_mumps_compiles AND NOT (_mumps_run STREQUAL "FAILED_TO_RUN")) break() endif() foreach(_pdep ${_mumps_potential_dependencies}) #message("CHECKING ${_pdep}") set(_add_pdep FALSE) if (NOT _mumps_compiles AND _out MATCHES "undefined reference.*${_mumps_dep_symbol_${_pdep}}") set(_add_pdep TRUE) #message("NEED COMPILE ${_pdep}") elseif(_mumps_run STREQUAL "FAILED_TO_RUN" AND DEFINED _mumps_run_dep_symbol_${_pdep} AND _run MATCHES "${_mumps_run_dep_symbol_${_pdep}}") set(_add_pdep TRUE) - #message("NEED RUN ${_pdep}") + #message("NEED RUN ${_pdep}") endif() if(_add_pdep) mumps_add_dependency(${_pdep} _libs ${_mumps_dep_comp_${_pdep}}) #message("ADDING ${_libs}") if(NOT _libs) message(FATAL_ERROR "MUMPS depends on ${_pdep} but no libraries where found") endif() list(APPEND _libraries_all ${_libs}) set(_retry_try_run TRUE) endif() endforeach() math(EXPR _retry_count "${_retry_count} + 1") endwhile() if(_retry_count GREATER 10) message(FATAL_ERROR "Do not know what to do to link with mumps on your system, I give up!") + message("Last compilation outputs: \n${_out} \n And last run output \n${_run}") endif() if(APPLE) # in doubt add some stuff because mumps was perhaps badly compiled mumps_add_dependency(pord _libs) list(APPEND _libraries_all ${_libs}) endif() set(MUMPS_LIBRARIES_ALL ${_libraries_all} PARENT_SCOPE) endfunction() mumps_find_dependencies() set(MUMPS_LIBRARIES ${MUMPS_LIBRARIES_ALL} CACHE INTERNAL "" FORCE) #=============================================================================== include(FindPackageHandleStandardArgs) if(CMAKE_VERSION VERSION_GREATER 2.8.12) if(MUMPS_INCLUDE_DIR) file(STRINGS ${MUMPS_INCLUDE_DIR}/dmumps_c.h _versions REGEX "^#define MUMPS_VERSION .*") foreach(_ver ${_versions}) string(REGEX MATCH "MUMPS_VERSION *\"([0-9.]+)\"" _tmp "${_ver}") set(_mumps_VERSION ${CMAKE_MATCH_1}) endforeach() set(MUMPS_VERSION "${_mumps_VERSION}" CACHE INTERNAL "") endif() find_package_handle_standard_args(Mumps REQUIRED_VARS ${_mumps_required_vars} MUMPS_INCLUDE_DIR VERSION_VAR MUMPS_VERSION ) else() find_package_handle_standard_args(Mumps DEFAULT_MSG ${_mumps_required_vars} MUMPS_INCLUDE_DIR) endif() diff --git a/examples/explicit/explicit_dynamic.cc b/examples/explicit/explicit_dynamic.cc index dbb57c5b1..39e61bb08 100644 --- a/examples/explicit/explicit_dynamic.cc +++ b/examples/explicit/explicit_dynamic.cc @@ -1,106 +1,107 @@ /** * @file explicit_dynamic.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * * @date creation: Sun Jul 12 2015 * @date last modification: Mon Jan 18 2016 * * @brief This code refers to the explicit dynamic example from the user manual * * * 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 "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize("material.dat", argc, argv); const UInt spatial_dimension = 3; const Real pulse_width = 2.; const Real A = 0.01; Real time_step; Real time_factor = 0.8; UInt max_steps = 1000; Mesh mesh(spatial_dimension); if (Communicator::getStaticCommunicator().whoAmI() == 0) mesh.read("bar.msh"); SolidMechanicsModel model(mesh); /// model initialization model.initFull(_analysis_method = _explicit_lumped_mass); time_step = model.getStableTimeStep(); std::cout << "Time Step = " << time_step * time_factor << "s (" << time_step << "s)" << std::endl; - model.setTimeStep(time_step * time_factor); + time_step *= time_factor; + model.setTimeStep(time_step); /// boundary and initial conditions Array & displacement = model.getDisplacement(); const Array & nodes = mesh.getNodes(); for (UInt n = 0; n < mesh.getNbNodes(); ++n) { Real x = nodes(n) - 2; // Sinus * Gaussian Real L = pulse_width; Real k = 0.1 * 2 * M_PI * 3 / L; displacement(n) = A * sin(k * x) * exp(-(k * x) * (k * x) / (L * L)); } std::ofstream energy; energy.open("energy.csv"); energy << "id,rtime,epot,ekin,tot" << std::endl; model.setBaseName("explicit_dynamic"); model.addDumpField("displacement"); model.addDumpField("velocity"); model.addDumpField("acceleration"); model.addDumpField("stress"); model.dump(); for (UInt s = 1; s <= max_steps; ++s) { model.solveStep(); Real epot = model.getEnergy("potential"); Real ekin = model.getEnergy("kinetic"); energy << s << "," << s * time_step << "," << epot << "," << ekin << "," << epot + ekin << "," << std::endl; if (s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; model.dump(); } energy.close(); finalize(); return EXIT_SUCCESS; } diff --git a/examples/python/CMakeLists.txt b/examples/python/CMakeLists.txt index 343e776e5..8f5758fd4 100644 --- a/examples/python/CMakeLists.txt +++ b/examples/python/CMakeLists.txt @@ -1,9 +1,10 @@ add_subdirectory(custom-material) add_subdirectory(dynamics) add_subdirectory(eigen_modes) add_subdirectory(plate-hole) add_subdirectory(stiffness_matrix) +add_subdirectory(structural_mechanics) package_add_files_to_package( examples/python/README.rst ) diff --git a/examples/python/structural_mechanics/CMakeLists.txt b/examples/python/structural_mechanics/CMakeLists.txt new file mode 100644 index 000000000..f8a6b2751 --- /dev/null +++ b/examples/python/structural_mechanics/CMakeLists.txt @@ -0,0 +1,14 @@ +register_example(structural_mechanics.py + SCRIPT structural_mechanics.py + PYTHON + ) + +register_example(structural_mechanics_softening.py + SCRIPT structural_mechanics.py + PYTHON + ) + +register_example(structural_mechanics_dynamics.py + SCRIPT structural_mechanics.py + PYTHON + ) diff --git a/examples/python/structural_mechanics/structural_mechanics.py b/examples/python/structural_mechanics/structural_mechanics.py new file mode 100644 index 000000000..8599a50a0 --- /dev/null +++ b/examples/python/structural_mechanics/structural_mechanics.py @@ -0,0 +1,135 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Test of Structural Mechanics +# In this example a beam, consisting of two elements, three nodes, is created. +# The left most node is fixed and a force is applied at the right most node. +import akantu as aka +import numpy +import matplotlib.pyplot as plt +import numpy as np + +# ### Creating the Mesh +# Create a mesh for the two dimensional case +beam = aka.Mesh(2) + +# We now create the connectivity array for the beam. +beam.addConnectivity(aka._bernoulli_beam_2) + +# We need a `MeshAccessor` in order to change the size of the mesh entities. +beamAcc = aka.MeshAccessor(beam) + +# Now we create the array to store the nodes and the connectivities and give them their size. +beamAcc.resizeConnectivity(2, aka._bernoulli_beam_2) +beamAcc.resizeNodes(3) + +Nodes = beam.getNodes() +Nodes[0, :] = [0., 0.] +Nodes[1, :] = [1., 0.] +Nodes[2, :] = [2., 0.] + +# #### Setting the Connections +Conn = beam.getConnectivity(aka._bernoulli_beam_2) +Conn[0, :] = [0, 1] +Conn[1, :] = [1, 2] + +# #### Ready +# We have to make the mesh ready. +beamAcc.makeReady() + + +# ### Creating the Model +model = aka.StructuralMechanicsModel(beam) + +# #### Setting up the Modell +# ##### Creating and Inserting the Materials +mat1 = aka.StructuralMaterial() +mat1.E = 1e9 +mat1.rho = 1. +mat1.I = 1. +mat1.Iz = 1. +mat1.Iy = 1. +mat1.A = 1. +mat1.GJ = 1. +model.addMaterial(mat1) + +mat2 = aka.StructuralMaterial() +mat2.E = 1e9 +mat2.rho = 1. +mat2.I = 1. +mat2.Iz = 1. +mat2.Iy = 1. +mat2.A = 1. +mat2.GJ = 1. +model.addMaterial(mat2) + +# ##### Initializing the Model +model.initFull(aka._implicit_dynamic) + +# ##### Assigning the Materials +materials = model.getElementMaterial(aka._bernoulli_beam_2) +materials[0][0] = 0 +materials[1][0] = 1 + +# ##### Setting Boundaries + +# Neumann +# Apply a force of `10` at the last (right most) node. +forces = model.getExternalForce() +forces[:] = 0 +forces[2, 0] = 100. + +# Dirichlets +# Block all dofs of the first node, since it is fixed. +# All other nodes have no restrictions +boundary = model.getBlockedDOFs() +boundary[0, :] = True +boundary[1, :] = False +boundary[2, :] = False + +# ### Solving the System + +# Set up the system +deltaT = 1e-10 +model.setTimeStep(deltaT) +solver = model.getNonLinearSolver() +solver.set("max_iterations", 100) +solver.set("threshold", 1e-8) +solver.set("convergence_type", aka.SolveConvergenceCriteria.solution) + +# Perform N time steps. +# At each step records the displacement of all three nodes in x direction. +N = 1000000 + +disp1 = np.zeros(N) +disp2 = np.zeros(N) +disp0 = np.zeros(N) +times = np.zeros(N) + +for i in range(N): + model.solveStep() + disp = model.getDisplacement() + disp0[i] = disp[0, 0] + disp1[i] = disp[1, 0] + disp2[i] = disp[2, 0] + times[i] = deltaT * i + +disps = [disp0, disp1, disp2] +maxMin = [-1.0, 1.0] + +for d in disps: + maxMin[0] = max(np.max(d), maxMin[0]) + maxMin[1] = min(np.min(d), maxMin[1]) + +plt.plot(disp1, times, color='g', label = "middle node") +plt.plot(disp2, times, color='b', label = "right node") + +plt.title("Displacement in $x$ of the nodes") +plt.ylabel("Time [S]") +plt.xlabel("displacement [m]") + +plt.xlim((maxMin[1] * 1.3, maxMin[0] * 1.1)) + +plt.legend() + +plt.show() diff --git a/examples/python/structural_mechanics/structural_mechanics_dynamics.py b/examples/python/structural_mechanics/structural_mechanics_dynamics.py new file mode 100644 index 000000000..b3d7cf850 --- /dev/null +++ b/examples/python/structural_mechanics/structural_mechanics_dynamics.py @@ -0,0 +1,166 @@ +#!/usr/bin/env python +# coding: utf-8 + +import akantu as aka +import numpy +import matplotlib.pyplot as plt +import numpy as np + +# ### Creating the Mesh +# Create a mesh for the two dimensional case +el_type = aka._bernoulli_beam_2 +beam = aka.Mesh(2) + +# We now create the connectivity array for the beam. +beam.addConnectivityType(el_type) + +# We need a `MeshAccessor` in order to change the size of the mesh entities. +beamAcc = aka.MeshAccessor(beam) + +# Now we create the array to store the nodes and the connectivities and give them their size. +nb_elem = 40 +L = 2 +beamAcc.resizeConnectivity(nb_elem, el_type) +beamAcc.resizeNodes(nb_elem + 1) + +# #### Setting the Nodes +Nodes = beam.getNodes() +length = L / nb_elem +Nodes[:, :] = 0. +Nodes[:, 0] = np.arange(nb_elem+1) * length + +# #### Setting the Connections +Conn = beam.getConnectivity(el_type) + +for e in range(nb_elem): + Conn[e, :] = [e, e + 1] + +# #### Ready +# We have to make the mesh ready. +beamAcc.makeReady() + +# ### Creating the Model +model = aka.StructuralMechanicsModel(beam) + +if el_type == aka._bernoulli_beam_3: + normal = beam.getDataReal("extra_normal", el_type) + + for e in range(nb_elem): + normal[e, :] = [0, 0, 1] + +# #### Setting up the Modell +# ##### Creating and Inserting the Materials +mat1 = aka.StructuralMaterial() +mat1.E = 1e9 +mat1.rho = 10. +mat1.I = 1. +mat1.Iz = 1. +mat1.Iy = 1. +mat1.A = 1. +mat1.GJ = 1. +model.addMaterial(mat1, 'mat1') + +# ##### Initializing the Model +model.initFull(aka.AnalysisMethod._implicit_dynamic) + +# ##### Assigning the Materials +materials = model.getElementMaterial(el_type) +materials[:, :] = 0 + +# ##### Setting Boundaries +# Neumann +F = 1e4 +no_print = int(nb_elem / 2) + +# Apply a force of `10` at the last (right most) node. +forces = model.getExternalForce() +forces[:, :] = 0 +forces[no_print, 1] = F + +# Dirichlets +# Block all dofs of the first node, since it is fixed. +# All other nodes have no restrictions +boundary = model.getBlockedDOFs() +boundary[:, :] = False + +boundary[0, 0] = True +boundary[0, 1] = True + +if el_type == aka._bernoulli_beam_3: + boundary[0, 2] = True + +boundary[nb_elem, 1] = True + +# ### Solving the System +# Set up the system +deltaT = 1e-6 +model.setTimeStep(deltaT) +solver = model.getNonLinearSolver() +solver.set("max_iterations", 100) +solver.set("threshold", 1e-8) +solver.set("convergence_type", aka.SolveConvergenceCriteria.solution) + +model.assembleMatrix("M") +M_ = model.getDOFManager().getMatrix("M") +M = aka.AkantuSparseMatrix(M_) + +model.assembleMatrix("K") +K_ = model.getDOFManager().getMatrix("K") +K = aka.AkantuSparseMatrix(K_) + +C_ = model.getDOFManager().getMatrix("C") +C_.add(M_, 0.00001) +C_.add(K_, 0.00001) + +def analytical_solution(time, L, rho, E, A, I, F): + omega = np.pi**2 / L**2 * np.sqrt(E * I / rho) + sum = 0. + N = 110 + for n in range(1, N, 2): + sum += (1. - np.cos(n * n * omega * time)) / n**4 + + return 2. * F * L**3 / np.pi**4 / E / I * sum + +# Perform N time steps. +# At each step records the displacement of all three nodes in x direction. +N = 900 + +mat1 = model.getMaterial('mat1') + +disp = model.getDisplacement() +velo = model.getVelocity() +disp[:, :] = 0. + +displs = np.zeros(N) + +ekin = np.zeros(N) +epot = np.zeros(N) +ework = np.zeros(N) +_ework = 0. + +for i in range(1, N): + model.solveStep() + displs[i] = disp[no_print, 1] + + _ework += F * velo[no_print, 1] * deltaT + + ekin[i] = model.getEnergy("kinetic") + epot[i] = model.getEnergy("potential") + ework[i] = _ework + + +def sol(x): + return analytical_solution(x, L, mat1.rho, mat1.E, + mat1.A, mat1.I, F) + + +times = np.arange(N) * deltaT +plt.plot(times, sol(times)) +plt.plot(times, displs) +plt.plot(times, displs - sol(times)) + +# What I do not fully understand is why the middle node first go backwards until it goes forward. +# I could imagine that there is some vibration, because everything is in rest. +np.max(displs - sol(times)) +plt.plot(times, ekin+epot) +plt.plot(times, ework) diff --git a/examples/python/structural_mechanics/structural_mechanics_softening.py b/examples/python/structural_mechanics/structural_mechanics_softening.py new file mode 100644 index 000000000..b196c8a71 --- /dev/null +++ b/examples/python/structural_mechanics/structural_mechanics_softening.py @@ -0,0 +1,188 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Test of Structural Mechanics +# In this test there is a beam consisting of three parts, all have the same materials. +# The left most node is fixed. +# On the right most node a force is applied in x direction. +# +# After a certain time, the material of the middle _element_ is waekened, lower Young's modulus. +# In each step the modulus is lowered by a coinstant factor. + +import akantu as aka +import numpy +import matplotlib.pyplot as plt +import numpy as np + + +# ### Creating the Mesh +# Create a mesh for the two dimensional case +beam = aka.Mesh(2) + +# We now create the connectivity array for the beam. +beam.addConnectivityType(aka._bernoulli_beam_2) + +# We need a `MeshAccessor` in order to change the size of the mesh entities. +beamAcc = aka.MeshAccessor(beam) + +# Now we create the array to store the nodes and the connectivities and give them their size. +beamAcc.resizeConnectivity(3, aka._bernoulli_beam_2) +beamAcc.resizeNodes(4) + +# #### Setting the Nodes +Nodes = beam.getNodes() +Nodes[0, :] = [0., 0.] +Nodes[1, :] = [1., 0.] +Nodes[2, :] = [2., 0.] +Nodes[3, :] = [3., 0.] + + +# #### Setting the Connections +Conn = beam.getConnectivity(aka._bernoulli_beam_2) +Conn[0, :] = [0, 1] +Conn[1, :] = [1, 2] +Conn[2, :] = [2, 3] + +#### Ready +# We have to make the mesh ready. +beamAcc.makeReady() + +# ### Creating the Model +model = aka.StructuralMechanicsModel(beam) + +# #### Setting up the Modell +# ##### Creating and Inserting the Materials +mat1 = aka.StructuralMaterial() +mat1.E = 1e9 +mat1.rho = 1. +mat1.I = 1. +mat1.Iz = 1. +mat1.Iy = 1. +mat1.A = 1. +mat1.GJ = 1. +mat1ID = model.addMaterial(mat1, 'mat1') + +mat2 = aka.StructuralMaterial() +mat2.E = 1e9 +mat2.rho = 1. +mat2.I = 1. +mat2.Iz = 1. +mat2.Iy = 1. +mat2.A = 1. +mat2.GJ = 1. +mat2ID = model.addMaterial(mat2, 'mat2') + +mat3 = aka.StructuralMaterial() +mat3.E = mat2.E / 100000 +mat3.rho = 1. +mat3.I = 1. +mat3.Iz = 1. +mat3.Iy = 1. +mat3.A = mat2.A / 100 +mat3.GJ = 1. +mat3ID = model.addMaterial(mat3, 'mat3') + +# ##### Initializing the Model +model.initFull(aka.AnalysisMethod._implicit_dynamic) + +# ##### Assigning the Materials +materials = model.getElementMaterial(aka._bernoulli_beam_2) + +materials[0][0] = mat1ID +materials[1][0] = mat2ID +materials[2][0] = mat1ID + + +# ##### Setting Boundaries +# Neumann +# Apply a force of `10` at the last (right most) node. +forces = model.getExternalForce() +forces[:] = 0 +forces[2, 0] = 100. + +# Dirichlets +# Block all dofs of the first node, since it is fixed. +# All other nodes have no restrictions +boundary = model.getBlockedDOFs() +boundary[0, :] = True +boundary[1, :] = False +boundary[2, :] = False +boundary[3, :] = False + +# ### Solving the System +# Set up the system +deltaT = 1e-9 +model.setTimeStep(deltaT) +solver = model.getNonLinearSolver() +solver.set("max_iterations", 100) +solver.set("threshold", 1e-8) +solver.set("convergence_type", aka.SolveConvergenceCriteria.solution) + +# Perform N time steps. +# At each step records the displacement of all three nodes in x direction. +N = 10000 * 60 + +disp0 = np.zeros(N) +disp1 = np.zeros(N) +disp2 = np.zeros(N) +disp3 = np.zeros(N) +times = np.zeros(N) +switchT = None +switchEnd = None + +softDuration = 1000 +SoftStart = (N // 2) - softDuration // 2 +SoftEnd = SoftStart + softDuration +if(softDuration > 0): + softFactor = (model.getMaterial('mat3ID').E / model.getMaterial('mat2ID').E) ** (1.0 / softDuration) + +mat2 = model.getMaterial('mat2') + +for i in range(N): + times[i] = deltaT * i + + if((SoftStart <= i <= SoftEnd) and (softDuration > 0)): + if switchT is None: + switchT = times[i] + elif(i == SoftEnd): + switchEnd = times[i] + # + mat2.E *= softFactor + # + + model.solveStep() + disp = model.getDisplacement() + disp0[i] = disp[0, 0] + disp1[i] = disp[1, 0] + disp2[i] = disp[2, 0] + disp3[i] = disp[3, 0] + +disps = [disp0, disp1, disp2, disp3] +maxMin = [-1.0, 1.0] + +for d in disps: + maxMin[0] = max(np.max(d), maxMin[0]) + maxMin[1] = min(np.min(d), maxMin[1]) + +#plt.plot(disp0, times, color='k', label = "left node (fix)") +plt.plot(disp1, times, color='g', label = "middle, left node") +plt.plot(disp2, times, color='g', linestyle = '--', label = "middle, right node") +plt.plot(disp3, times, color='b', label = "right node") + +if(softDuration > 0): + plt.plot((maxMin[1], maxMin[0]), (switchT, switchT),) + plt.plot((maxMin[1], maxMin[0]), (switchEnd, switchEnd), ) + +plt.title("Displacement in $x$ of the nodes") +plt.ylabel("Time [S]") +plt.xlabel("displacement [m]") +plt.xlim((maxMin[1] * 1.3, maxMin[0] * 1.1)) +plt.legend() +plt.show() + +# If the softening is disabled, then the displacement looks wierd. +# Because the displacement first increases and then decreases. +# In this case `softDuration > 0` holds. +# +# However if the softening is enabled, it looks rather good. +# The left middle node will start to vibrate, because it is not pulled in the other direction. diff --git a/examples/structural_mechanics/CMakeLists.txt b/examples/structural_mechanics/CMakeLists.txt index 75d257f84..fb1688a8d 100644 --- a/examples/structural_mechanics/CMakeLists.txt +++ b/examples/structural_mechanics/CMakeLists.txt @@ -1,37 +1,37 @@ #=============================================================================== # @file CMakeLists.txt # # @author Fabian Barras # # @date creation: Mon Jan 18 2016 # @date last modification: Tue Jan 19 2016 # # @brief configuration for structural mechanics example # # @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 . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== -#register_example(bernoulli_beam_2_exemple -# SOURCES bernoulli_beam_2_exemple.cc -# ) +register_example(bernoulli_beam_2_example + SOURCES bernoulli_beam_2_example.cc + ) diff --git a/examples/structural_mechanics/bernoulli_beam_2_exemple.cc b/examples/structural_mechanics/bernoulli_beam_2_example.cc similarity index 78% rename from examples/structural_mechanics/bernoulli_beam_2_exemple.cc rename to examples/structural_mechanics/bernoulli_beam_2_example.cc index b825f93e7..c2f250518 100644 --- a/examples/structural_mechanics/bernoulli_beam_2_exemple.cc +++ b/examples/structural_mechanics/bernoulli_beam_2_example.cc @@ -1,164 +1,151 @@ /** * @file bernoulli_beam_2_exemple.cc * * @author Fabian Barras * * @date creation: Mon Jan 18 2016 * * @brief Computation of the analytical exemple 1.1 in the TGC vol 6 * * * 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 "structural_mechanics_model.hh" +#include "mesh_accessor.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #define TYPE _bernoulli_beam_2 using namespace akantu; // Linear load function -static void lin_load(double * position, double * load, - __attribute__((unused)) Real * normal, - __attribute__((unused)) UInt surface_id) { - memset(load, 0, sizeof(Real) * 3); - if (position[0] <= 10) { - load[1] = -6000; +static void lin_load(const Array & nodes, Array& load) { + for(auto &&data : zip(make_view(nodes, 2), make_view(load, 2))) { + if (std::get<0>(data)[_y] <= 10) { + std::get<1>(data)[_y] = -6000; + } } } /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); // Defining the mesh Mesh beams(2); UInt nb_nodes = 3; UInt nb_nodes_1 = 1; UInt nb_nodes_2 = nb_nodes - nb_nodes_1 - 1; UInt nb_element = nb_nodes - 1; - Array & nodes = const_cast &>(beams.getNodes()); + MeshAccessor mesh_accessor(beams); + Array & nodes = mesh_accessor.getNodes(); nodes.resize(nb_nodes); beams.addConnectivityType(_bernoulli_beam_2); - Array & connectivity = - const_cast &>(beams.getConnectivity(_bernoulli_beam_2)); + Array & connectivity = mesh_accessor.getConnectivity(_bernoulli_beam_2); connectivity.resize(nb_element); for (UInt i = 0; i < nb_nodes; ++i) { nodes(i, 1) = 0; } for (UInt i = 0; i < nb_nodes_1; ++i) { nodes(i, 0) = 10. * i / ((Real)nb_nodes_1); } nodes(nb_nodes_1, 0) = 10; for (UInt i = 0; i < nb_nodes_2; ++i) { nodes(nb_nodes_1 + i + 1, 0) = 10 + 8. * (i + 1) / ((Real)nb_nodes_2); } for (UInt i = 0; i < nb_element; ++i) { connectivity(i, 0) = i; connectivity(i, 1) = i + 1; } // Defining the materials StructuralMechanicsModel model(beams); StructuralMaterial mat1; mat1.E = 3e10; mat1.I = 0.0025; mat1.A = 0.01; model.addMaterial(mat1); StructuralMaterial mat2; mat2.E = 3e10; mat2.I = 0.00128; mat2.A = 0.01; model.addMaterial(mat2); // Defining the forces model.initFull(); const Real M = -3600; // Momentum at 3 - Array & forces = model.getForce(); + Array & forces = model.getExternalForce(); Array & displacement = model.getDisplacement(); Array & boundary = model.getBlockedDOFs(); const Array & N_M = model.getStress(_bernoulli_beam_2); Array & element_material = model.getElementMaterial(_bernoulli_beam_2); forces.zero(); displacement.zero(); for (UInt i = 0; i < nb_nodes_2; ++i) { element_material(i + nb_nodes_1) = 1; } forces(nb_nodes - 1, 2) += M; - model.computeForcesFromFunction<_bernoulli_beam_2>(lin_load, _bft_traction); + Array load(nodes.size(), 2); + lin_load(nodes, load); + + model.computeForcesByGlobalTractionArray<_bernoulli_beam_2>(load); // Defining the boundary conditions boundary(0, 0) = true; boundary(0, 1) = true; boundary(0, 2) = true; boundary(nb_nodes_1, 1) = true; boundary(nb_nodes - 1, 1) = true; - // Solve - Real error; - model.assembleStiffnessMatrix(); - - UInt count = 0; model.addDumpFieldVector("displacement"); model.addDumpField("rotation"); model.addDumpFieldVector("force"); model.addDumpField("momentum"); - do { - if (count != 0) - std::cerr << count << " - " << error << std::endl; - model.updateResidual(); - model.solve(); - count++; - } while (!model.testConvergenceIncrement(1e-10, error) && count < 10); - std::cerr << count << " - " << error << std::endl; - - /* -------------------------------------------------------------------------- - */ - // Post-Processing + model.solveStep(); - model.computeStresses(); + // Post-Processing std::cout << " d1 = " << displacement(nb_nodes_1, 2) << std::endl; std::cout << " d2 = " << displacement(nb_nodes - 1, 2) << std::endl; std::cout << " M1 = " << N_M(0, 1) << std::endl; std::cout << " M2 = " << N_M(2 * (nb_nodes - 2), 1) << std::endl; model.dump(); finalize(); } diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 5ac5d0380..76a710dc2 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -1,109 +1,118 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date creation: Fri Dec 12 2014 # @date last modification: Mon Jan 18 2016 # # @brief CMake file for the python wrapping of akantu # # @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 . # #=============================================================================== package_get_all_include_directories( AKANTU_LIBRARY_INCLUDE_DIRS ) package_get_all_external_informations( PRIVATE_INCLUDE AKANTU_PRIVATE_EXTERNAL_INCLUDE_DIR INTERFACE_INCLUDE AKANTU_INTERFACE_EXTERNAL_INCLUDE_DIR LIBRARIES AKANTU_EXTERNAL_LIBRARIES ) set(PYAKANTU_SRCS py_aka_common.cc py_aka_error.cc py_akantu.cc py_boundary_conditions.cc py_fe_engine.cc py_group_manager.cc py_mesh.cc py_model.cc py_parser.cc + py_solver.cc ) package_is_activated(iohelper _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_dumpable.cc ) endif() package_is_activated(solid_mechanics _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_solid_mechanics_model.cc py_material.cc ) endif() package_is_activated(cohesive_element _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_solid_mechanics_model_cohesive.cc ) endif() package_is_activated(heat_transfer _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_heat_transfer_model.cc ) endif() +package_is_activated(structural_mechanics _is_activated) +if (_is_activated) + list(APPEND PYAKANTU_SRCS + py_structural_mechanics_model.cc + ) +endif() + + pybind11_add_module(py11_akantu ${PYAKANTU_SRCS}) # to avoid compilation warnings from pybind11 target_include_directories(py11_akantu SYSTEM BEFORE PRIVATE ${PYBIND11_INCLUDE_DIR} PRIVATE ${pybind11_INCLUDE_DIR} PRIVATE ${PYTHON_INCLUDE_DIRS}) target_link_libraries(py11_akantu PUBLIC akantu) set_target_properties(py11_akantu PROPERTIES DEBUG_POSTFIX "") set(_python_install_dir ${CMAKE_INSTALL_LIBDIR}/python${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}/site-packages) install(TARGETS py11_akantu LIBRARY DESTINATION ${_python_install_dir}) install(DIRECTORY akantu DESTINATION ${_python_install_dir} FILES_MATCHING PATTERN "*.py") configure_file(setup.py.in "${CMAKE_CURRENT_BINARY_DIR}/setup.py" @ONLY) diff --git a/python/py_akantu.cc b/python/py_akantu.cc index 666a20e7b..d3d3a93a8 100644 --- a/python/py_akantu.cc +++ b/python/py_akantu.cc @@ -1,103 +1,113 @@ /* -------------------------------------------------------------------------- */ #include "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "py_aka_common.hh" #include "py_aka_error.hh" #include "py_boundary_conditions.hh" #include "py_fe_engine.hh" #include "py_group_manager.hh" #include "py_mesh.hh" #include "py_model.hh" #include "py_parser.hh" +#include "py_solver.hh" #if defined(AKANTU_USE_IOHELPER) #include "py_dumpable.hh" #endif #if defined(AKANTU_SOLID_MECHANICS) #include "py_material.hh" #include "py_solid_mechanics_model.hh" #endif #if defined(AKANTU_HEAT_TRANSFER) #include "py_heat_transfer_model.hh" #endif #if defined(AKANTU_COHESIVE_ELEMENT) #include "py_solid_mechanics_model_cohesive.hh" #endif + +#if defined(AKANTU_STRUCTURAL_MECHANICS) +#include "py_structural_mechanics_model.hh" +#endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; namespace akantu { void register_all(pybind11::module & mod) { register_initialize(mod); register_enums(mod); register_error(mod); register_functions(mod); register_parser(mod); + register_solvers(mod); register_group_manager(mod); #if defined(AKANTU_USE_IOHELPER) register_dumpable(mod); #endif register_mesh(mod); register_fe_engine(mod); register_boundary_conditions(mod); register_model(mod); #if defined(AKANTU_HEAT_TRANSFER) register_heat_transfer_model(mod); #endif #if defined(AKANTU_SOLID_MECHANICS) register_solid_mechanics_model(mod); register_material(mod); #endif #if defined(AKANTU_COHESIVE_ELEMENT) register_solid_mechanics_model_cohesive(mod); #endif + +#if defined(AKANTU_STRUCTURAL_MECHANICS) + register_structural_mechanics_model(mod); +#endif } } // namespace akantu /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ PYBIND11_MODULE(py11_akantu, mod) { mod.doc() = "Akantu python interface"; static py::exception akantu_exception(mod, "Exception"); py::register_exception_translator([](std::exception_ptr ptr) { try { if (ptr) { std::rethrow_exception(ptr); } } catch (akantu::debug::Exception & e) { if (akantu::debug::debugger.printBacktrace()) { akantu::debug::printBacktrace(); } akantu_exception(e.info().c_str()); } }); akantu::register_all(mod); mod.def("has_mpi", []() { #if defined(AKANTU_USE_MPI) return true; #else return false; #endif }); } // Module akantu diff --git a/python/py_material.cc b/python/py_material.cc index cf0c02f76..a35063a47 100644 --- a/python/py_material.cc +++ b/python/py_material.cc @@ -1,231 +1,207 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #if defined(AKANTU_COHESIVE_ELEMENT) #include #endif /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ #if not defined(PYBIND11_OVERRIDE) #define PYBIND11_OVERRIDE PYBIND11_OVERLOAD #define PYBIND11_OVERRIDE_PURE PYBIND11_OVERLOAD_PURE #endif namespace akantu { template class PyMaterial : public _Material { public: /* Inherit the constructors */ using _Material::_Material; ~PyMaterial() override = default; void initMaterial() override { PYBIND11_OVERRIDE(void, _Material, initMaterial, ); // NOLINT }; void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override { PYBIND11_OVERRIDE_PURE(void, _Material, computeStress, el_type, ghost_type); } void computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override { PYBIND11_OVERRIDE(void, _Material, computeTangentModuli, el_type, tangent_matrix, ghost_type); } void computePotentialEnergy(ElementType el_type) override { PYBIND11_OVERRIDE(void, _Material, computePotentialEnergy, el_type); } Real getPushWaveSpeed(const Element & element) const override { PYBIND11_OVERRIDE(Real, _Material, getPushWaveSpeed, element); } Real getShearWaveSpeed(const Element & element) const override { PYBIND11_OVERRIDE(Real, _Material, getShearWaveSpeed, element); } template void registerInternal(const std::string & name, UInt nb_component) { auto && internal = std::make_shared>(name, *this); AKANTU_DEBUG_INFO("alloc internal " << name << " " << &this->internals[name]); internal->initialize(nb_component); this->internals[name] = internal; } protected: std::map> internals; }; /* -------------------------------------------------------------------------- */ template -void register_element_type_map_array(py::module & mod, +void register_internal_field(py::module & mod, const std::string & name) { - py::class_, std::shared_ptr>>( - mod, ("ElementTypeMapArray" + name).c_str()) - .def( - "__call__", - [](ElementTypeMapArray & self, ElementType type, - GhostType ghost_type) -> decltype(auto) { - return self(type, ghost_type); - }, - py::arg("type"), py::arg("ghost_type") = _not_ghost, - py::return_value_policy::reference) - .def( - "elementTypes", - [](ElementTypeMapArray & self, UInt _dim, GhostType _ghost_type, - ElementKind _kind) -> decltype(auto) { - auto types = self.elementTypes(_dim, _ghost_type, _kind); - std::vector _types; - for (auto && t : types) { - _types.push_back(t); - } - return _types; - }, - py::arg("dim") = _all_dimensions, py::arg("ghost_type") = _not_ghost, - py::arg("kind") = _ek_regular); - py::class_, ElementTypeMapArray, std::shared_ptr>>( mod, ("InternalField" + name).c_str()); } /* -------------------------------------------------------------------------- */ template void define_material(py::module & mod, const std::string & name) { py::class_<_Material, PyMaterial<_Material>, Parsable>( mod, name.c_str(), py::multiple_inheritance()) .def(py::init()) .def( "getGradU", [](Material & self, ElementType el_type, GhostType ghost_type = _not_ghost) -> decltype(auto) { return self.getGradU(el_type, ghost_type); }, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getStress", [](Material & self, ElementType el_type, GhostType ghost_type = _not_ghost) -> decltype(auto) { return self.getStress(el_type, ghost_type); }, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getPotentialEnergy", [](Material & self, ElementType el_type) -> decltype(auto) { return self.getPotentialEnergy(el_type); }, py::return_value_policy::reference) .def("initMaterial", &Material::initMaterial) .def("getModel", &Material::getModel) .def("registerInternalReal", [](Material & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .registerInternal(name, nb_component); }) .def("registerInternalUInt", [](Material & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .registerInternal(name, nb_component); }) .def( "getInternalReal", [](Material & self, const ID & id) -> decltype(auto) { return self.getInternal(id); }, py::arg("id"), py::return_value_policy::reference) .def( "getInternalUInt", [](Material & self, const ID & id) -> decltype(auto) { return self.getInternal(id); }, py::arg("id"), py::return_value_policy::reference) .def( "getElementFilter", [](Material & self) -> decltype(auto) { return self.getElementFilter(); }, py::return_value_policy::reference); } /* -------------------------------------------------------------------------- */ void register_material(py::module & mod) { py::class_(mod, "MaterialFactory") .def_static( "getInstance", []() -> MaterialFactory & { return Material::getFactory(); }, py::return_value_policy::reference) .def("registerAllocator", [](MaterialFactory & self, const std::string id, py::function func) { self.registerAllocator( id, [func, id](UInt dim, const ID & /*unused*/, SolidMechanicsModel & model, const ID & option) -> std::unique_ptr { py::object obj = func(dim, id, model, option); auto & ptr = py::cast(obj); obj.release(); return std::unique_ptr(&ptr); }); }); - register_element_type_map_array(mod, "Real"); - register_element_type_map_array(mod, "UInt"); + register_internal_field(mod, "Real"); + register_internal_field(mod, "UInt"); define_material(mod, "Material"); py::class_>( mod, "MaterialSelector") .def("setFallback", [](MaterialSelector & self, UInt f) { self.setFallback(f); }) .def("setFallback", [](MaterialSelector & self, const std::shared_ptr & fallback_selector) { self.setFallback(fallback_selector); }) .def("setFallback", [](MaterialSelector & self, MaterialSelector & fallback_selector) { self.setFallback(fallback_selector); }); py::class_, MaterialSelector, std::shared_ptr>>( mod, "MeshDataMaterialSelectorString") .def(py::init(), py::arg("name"), py::arg("model"), py::arg("first_index") = 1); #if defined(AKANTU_COHESIVE_ELEMENT) py::class_>( mod, "DefaultMaterialCohesiveSelector") .def(py::init()); py::class_>( mod, "MeshDataMaterialCohesiveSelector") .def(py::init()); py::class_>( mod, "MaterialCohesiveRulesSelector") .def(py::init(), py::arg("model"), py::arg("rules"), py::arg("mesh_data_id") = "physical_names"); #endif } } // namespace akantu diff --git a/python/py_mesh.cc b/python/py_mesh.cc index 928887e23..c67e422e5 100644 --- a/python/py_mesh.cc +++ b/python/py_mesh.cc @@ -1,71 +1,147 @@ /* -------------------------------------------------------------------------- */ #include "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include +#include #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { +/* -------------------------------------------------------------------------- */ +template +void register_element_type_map_array(py::module & mod, + const std::string & name) { + py::class_, std::shared_ptr>>( + mod, ("ElementTypeMapArray" + name).c_str()) + .def( + "__call__", + [](ElementTypeMapArray & self, ElementType type, + GhostType ghost_type) -> decltype(auto) { + return self(type, ghost_type); + }, + py::arg("type"), py::arg("ghost_type") = _not_ghost, + py::return_value_policy::reference) + .def( + "elementTypes", + [](ElementTypeMapArray & self, UInt _dim, GhostType _ghost_type, + ElementKind _kind) -> decltype(auto) { + auto types = self.elementTypes(_dim, _ghost_type, _kind); + std::vector _types; + for (auto && t : types) { + _types.push_back(t); + } + return _types; + }, + py::arg("dim") = _all_dimensions, py::arg("ghost_type") = _not_ghost, + py::arg("kind") = _ek_regular); +} + /* -------------------------------------------------------------------------- */ void register_mesh(py::module & mod) { + + register_element_type_map_array(mod, "Real"); + register_element_type_map_array(mod, "UInt"); + py::class_(mod, "MeshData") .def( "getElementalDataUInt", - [](MeshData & _this, const ID & name) -> ElementTypeMapArray & { + [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, + py::return_value_policy::reference) + .def( + "getElementalDataReal", + [](MeshData & _this, const ID & name) -> decltype(auto) { + return _this.getElementalData(name); + }, py::return_value_policy::reference); py::class_(mod, "Mesh", py::multiple_inheritance()) .def(py::init(), py::arg("spatial_dimension"), py::arg("id") = "mesh") .def("read", &Mesh::read, py::arg("filename"), py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file") .def( "getNodes", [](Mesh & self) -> decltype(auto) { return self.getNodes(); }, py::return_value_policy::reference) .def("getNbNodes", &Mesh::getNbNodes) .def( "getConnectivity", [](Mesh & self, ElementType type) -> decltype(auto) { return self.getConnectivity(type); }, py::return_value_policy::reference) + .def( + "addConnectivityType", + [](Mesh & self, ElementType type, GhostType ghost_type) -> void { + self.addConnectivityType(type, ghost_type); + }, + py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("distribute", [](Mesh & self) { self.distribute(); }) .def("makePeriodic", [](Mesh & self, const SpatialDirection & direction) { self.makePeriodic(direction); }) .def( "getNbElement", - [](Mesh & self, const UInt spatial_dimension, - GhostType ghost_type, ElementKind kind) { + [](Mesh & self, const UInt spatial_dimension, GhostType ghost_type, + ElementKind kind) { return self.getNbElement(spatial_dimension, ghost_type, kind); }, py::arg("spatial_dimension") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_not_defined) .def( "getNbElement", - [](Mesh & self, ElementType type, - GhostType ghost_type) { + [](Mesh & self, ElementType type, GhostType ghost_type) { return self.getNbElement(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) - .def_static("getSpatialDimension", [](ElementType & type) { - return Mesh::getSpatialDimension(type); - }); + .def_static( + "getSpatialDimension", + [](ElementType & type) { return Mesh::getSpatialDimension(type); }) + .def( + "getDataReal", + [](Mesh & _this, const ID & name, ElementType type, + GhostType ghost_type) -> decltype(auto) { + return _this.getData(name, type, ghost_type); + }, + py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost, + py::return_value_policy::reference) + .def( + "hasDataReal", + [](Mesh & _this, const ID & name, ElementType type, + GhostType ghost_type) -> bool { + return _this.hasData(name, type, ghost_type); + }, + py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost); /* ------------------------------------------------------------------------ */ py::class_(mod, "MeshUtils") .def_static("buildFacets", &MeshUtils::buildFacets); + + py::class_(mod, "MeshAccessor") + .def(py::init(), py::arg("mesh")) + .def( + "resizeConnectivity", + [](MeshAccessor & self, UInt new_size, ElementType type, GhostType gt) + -> void { self.resizeConnectivity(new_size, type, gt); }, + py::arg("new_size"), py::arg("type"), + py::arg("ghost_type") = _not_ghost) + .def( + "resizeNodes", + [](MeshAccessor & self, UInt new_size) -> void { + self.resizeNodes(new_size); + }, + py::arg("new_size")) + .def("makeReady", &MeshAccessor::makeReady); } } // namespace akantu diff --git a/python/py_model.cc b/python/py_model.cc index 32682d8a3..795b30ac9 100644 --- a/python/py_model.cc +++ b/python/py_model.cc @@ -1,80 +1,89 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_model(py::module & mod) { py::class_(mod, "DOFManager") .def("getMatrix", &DOFManager::getMatrix, py::return_value_policy::reference) .def( "getNewMatrix", [](DOFManager & self, const std::string & name, const std::string & matrix_to_copy_id) -> decltype(auto) { return self.getNewMatrix(name, matrix_to_copy_id); }, py::return_value_policy::reference) .def( "getResidual", [](DOFManager & self) -> decltype(auto) { return self.getResidual(); }, py::return_value_policy::reference) .def("getArrayPerDOFs", &DOFManager::getArrayPerDOFs) + .def( + "hasMatrix", + [](DOFManager & self, const ID & name) -> bool { + return self.hasMatrix(name); + }, + py::arg("name")) .def("assembleToResidual", &DOFManager::assembleToResidual); py::class_(mod, "NonLinearSolver") .def( "set", [](NonLinearSolver & self, const std::string & id, const Real & val) { if (id == "max_iterations") { self.set(id, int(val)); } else { self.set(id, val); } }) .def("set", [](NonLinearSolver & self, const std::string & id, const SolveConvergenceCriteria & val) { self.set(id, val); }); py::class_(mod, "ModelSolver", py::multiple_inheritance()) .def("getNonLinearSolver", (NonLinearSolver & (ModelSolver::*)(const ID &)) & ModelSolver::getNonLinearSolver, py::arg("solver_id") = "", py::return_value_policy::reference) .def("solveStep", [](ModelSolver & self) { self.solveStep(); }) .def("solveStep", [](ModelSolver & self, const ID & solver_id) { self.solveStep(solver_id); }); py::class_(mod, "Model", py::multiple_inheritance()) .def("setBaseName", &Model::setBaseName) .def("getFEEngine", &Model::getFEEngine, py::arg("name") = "", py::return_value_policy::reference) .def("getFEEngineBoundary", &Model::getFEEngine, py::arg("name") = "", py::return_value_policy::reference) .def("addDumpFieldVector", &Model::addDumpFieldVector) .def("addDumpField", &Model::addDumpField) .def("setBaseNameToDumper", &Model::setBaseNameToDumper) .def("addDumpFieldVectorToDumper", &Model::addDumpFieldVectorToDumper) .def("addDumpFieldToDumper", &Model::addDumpFieldToDumper) - .def("dump", &Model::dump) + .def("dump", py::overload_cast<>(&Model::dump)) + .def("dump", py::overload_cast(&Model::dump)) + .def("dump", py::overload_cast(&Model::dump)) .def("initNewSolver", &Model::initNewSolver) .def("getDOFManager", &Model::getDOFManager, - py::return_value_policy::reference); + py::return_value_policy::reference) + .def("assembleMatrix", &Model::assembleMatrix); } } // namespace akantu diff --git a/python/py_solver.cc b/python/py_solver.cc index 6e292867d..a0146642d 100644 --- a/python/py_solver.cc +++ b/python/py_solver.cc @@ -1,45 +1,53 @@ /* -------------------------------------------------------------------------- */ +#include "py_solver.hh" #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ -#include "py_solver.h" -/* -------------------------------------------------------------------------- */#include +#include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -void register_solver(py::module & mod) { +void register_solvers(py::module & mod) { py::class_(mod, "SparseMatrix") .def("getMatrixType", &SparseMatrix::getMatrixType) .def("size", &SparseMatrix::size) .def("zero", &SparseMatrix::zero) .def("saveProfile", &SparseMatrix::saveProfile) .def("saveMatrix", &SparseMatrix::saveMatrix) .def( "add", [](SparseMatrix & self, UInt i, UInt j) { self.add(i, j); }, "Add entry in the profile") .def( "add", [](SparseMatrix & self, UInt i, UInt j, Real value) { self.add(i, j, value); }, "Add the value to the matrix") + .def( + "add", + [](SparseMatrix & self, SparseMatrix & A, Real alpha) { + self.add(A, alpha); + }, + "Add a matrix to the matrix", py::arg("A"), py::arg("alpha") = 1.) .def("__call__", [](const SparseMatrix & self, UInt i, UInt j) { return self(i, j); }); py::class_(mod, "SparseMatrixAIJ") .def("getIRN", &SparseMatrixAIJ::getIRN) .def("getJCN", &SparseMatrixAIJ::getJCN) .def("getA", &SparseMatrixAIJ::getA); py::class_(mod, "SolverVector"); } + +} // namespace akantu diff --git a/python/py_solver.hh b/python/py_solver.hh index 892606998..41e1124e8 100644 --- a/python/py_solver.hh +++ b/python/py_solver.hh @@ -1,12 +1,12 @@ #include #ifndef AKANTU_PY_AKA_SOLVER_HH_ #define AKANTU_PY_AKA_SOLVER_HH_ namespace akantu { -void register_solver(pybind11::module & mod); +void register_solvers(pybind11::module & mod); } #endif diff --git a/python/py_structural_mechanics_model.cc b/python/py_structural_mechanics_model.cc new file mode 100644 index 000000000..4e3faa765 --- /dev/null +++ b/python/py_structural_mechanics_model.cc @@ -0,0 +1,130 @@ +/* -------------------------------------------------------------------------- */ +#include "py_aka_array.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ +namespace py = pybind11; +/* -------------------------------------------------------------------------- */ + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +#define def_deprecated(func_name, mesg) \ + def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) + +#define def_function_nocopy(func_name) \ + def( \ + #func_name, \ + [](StructuralMechanicsModel & self) -> decltype(auto) { \ + return self.func_name(); \ + }, \ + py::return_value_policy::reference) + +#define def_function_(func_name) \ + def(#func_name, [](StructuralMechanicsModel & self) -> decltype(auto) { \ + return self.func_name(); \ + }) + +#define def_plainmember(M) def_readwrite(#M, &StructuralMaterial::M) +/* -------------------------------------------------------------------------- */ + +void register_structural_mechanics_model(pybind11::module & mod) { + /* First we have to register the material class + * The wrapper aims to mimic the behaviour of the real material. + */ + py::class_(mod, "StructuralMaterial") + .def(py::init<>()) + .def(py::init()) + .def_plainmember(E) + .def_plainmember(A) + .def_plainmember(I) + .def_plainmember(Iz) + .def_plainmember(Iy) + .def_plainmember(GJ) + .def_plainmember(rho) + .def_plainmember(t) + .def_plainmember(nu); + + /* Now we create the structural model wrapper + * Note that this is basically a port from the solid mechanic part. + */ + py::class_(mod, "StructuralMechanicsModel") + .def(py::init(), + py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, + py::arg("id") = "structural_mechanics_model", + py::arg("memory_id") = 0) + .def( + "initFull", + [](StructuralMechanicsModel & self, + const AnalysisMethod & analysis_method) -> void { + self.initFull(_analysis_method = analysis_method); + }, + py::arg("_analysis_method")) + .def("initFull", + [](StructuralMechanicsModel & self) -> void { self.initFull(); }) + .def_function_nocopy(getExternalForce) + .def_function_nocopy(getDisplacement) + .def_function_nocopy(getInternalForce) + .def_function_nocopy(getVelocity) + .def_function_nocopy(getAcceleration) + .def_function_nocopy(getInternalForce) + .def_function_nocopy(getBlockedDOFs) + .def_function_nocopy(getMesh) + + .def("setTimeStep", &StructuralMechanicsModel::setTimeStep, + py::arg("time_step"), py::arg("solver_id") = "") + .def( + "getElementMaterial", + [](StructuralMechanicsModel & self, const ElementType & type, + GhostType ghost_type) -> decltype(auto) { + return self.getElementMaterial(type, ghost_type); + }, + "This function returns the map that maps elements to materials.", + py::arg("type"), py::arg("ghost_type") = _not_ghost, + py::return_value_policy::reference) + .def( + "getMaterialByElement", + [](StructuralMechanicsModel & self, Element element) + -> decltype(auto) { return self.getMaterialByElement(element); }, + "This function returns the `StructuralMaterial` instance that is " + "associated with element `element`.", + py::arg("element"), py::return_value_policy::reference) + .def( + "addMaterial", + [](StructuralMechanicsModel & self, StructuralMaterial & mat, + const ID & name) -> UInt { return self.addMaterial(mat, name); }, + "This function adds the `StructuralMaterial` `mat` to `self`." + " The function returns the ID of the new material.", + py::arg("mat"), py::arg("name") = "") + .def( + "getMaterial", + [](StructuralMechanicsModel & self, UInt material_index) + -> decltype(auto) { return self.getMaterial(material_index); }, + "This function returns the `i`th material of `self`." + " Note a reference is returned which allows to modify the material inside `self`.", + py::arg("i"), + py::return_value_policy::reference) + .def( + "getMaterial", + [](StructuralMechanicsModel & self, const ID & name) + -> decltype(auto) { return self.getMaterial(name); }, + "This function returns the material with name `i` of `self`." + " Note a reference is returned which allows to modify the material inside `self`.", + py::arg("i"), + py::return_value_policy::reference) + .def( + "getNbMaterials", + [](StructuralMechanicsModel & self) { return self.getNbMaterials(); }, + "Returns the number of different materials inside `self`.") + .def("getKineticEnergy", &StructuralMechanicsModel::getKineticEnergy, + "Compute kinetic energy") + .def("getPotentialEnergy", &StructuralMechanicsModel::getPotentialEnergy, + "Compute potential energy") + .def("getEnergy", &StructuralMechanicsModel::getEnergy, + "Compute the specified energy"); + +} // End: register structural mechanical model + +} // namespace akantu diff --git a/python/py_structural_mechanics_model.hh b/python/py_structural_mechanics_model.hh new file mode 100644 index 000000000..a859cd1f4 --- /dev/null +++ b/python/py_structural_mechanics_model.hh @@ -0,0 +1,12 @@ +#include + +#ifndef AKANTU_PY_STRUCTURAL_MECHANICS_MODEL_HH_ +#define AKANTU_PY_STRUCTURAL_MECHANICS_MODEL_HH_ + +namespace akantu { + +void register_structural_mechanics_model(pybind11::module & mod); + +} // namespace akantu + +#endif // AKANTU_PY_STRUCTURAL_MECHANICS_MODEL_HH_ diff --git a/src/fe_engine/element_class_structural.hh b/src/fe_engine/element_class_structural.hh index 2d9ae92a9..76732da6a 100644 --- a/src/fe_engine/element_class_structural.hh +++ b/src/fe_engine/element_class_structural.hh @@ -1,253 +1,274 @@ /** * @file element_class_structural.hh * * @author Fabian Barras * @author Lucas Frerot * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Thu Feb 21 2013 * @date last modification: Tue Feb 20 2018 * * @brief Specialization of the element classes for structural elements * * * Copyright (©) 2014-2018 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 "element_class.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_CLASS_STRUCTURAL_HH_ #define AKANTU_ELEMENT_CLASS_STRUCTURAL_HH_ namespace akantu { /// Macro to generate the InterpolationProperty structures for different /// interpolation types #define AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY( \ itp_type, itp_geom_type, ndof, nb_stress, nb_dnds_cols) \ template <> struct InterpolationProperty { \ static const InterpolationKind kind{_itk_structural}; \ static const UInt nb_nodes_per_element{ \ InterpolationProperty::nb_nodes_per_element}; \ static const InterpolationType itp_geometry_type{itp_geom_type}; \ static const UInt natural_space_dimension{ \ InterpolationProperty::natural_space_dimension}; \ static const UInt nb_degree_of_freedom{ndof}; \ static const UInt nb_stress_components{nb_stress}; \ static const UInt dnds_columns{nb_dnds_cols}; \ } /* -------------------------------------------------------------------------- */ template class InterpolationElement { public: using interpolation_property = InterpolationProperty; /// compute the shape values for a given set of points in natural coordinates static inline void computeShapes(const Matrix & natural_coord, const Matrix & real_coord, - Tensor3 & N) { + const Matrix & T, Tensor3 & Ns) { for (UInt i = 0; i < natural_coord.cols(); ++i) { - Matrix n_t = N(i); - computeShapes(natural_coord(i), real_coord, n_t); + Matrix N_T = Ns(i); + Matrix N(N_T.rows(), N_T.cols()); + computeShapes(natural_coord(i), real_coord, N); + N_T.mul(N, T); } } /// compute the shape values for a given point in natural coordinates static inline void computeShapes(const Vector & natural_coord, const Matrix & real_coord, Matrix & N); + static inline void computeShapesMass(const Matrix & natural_coords, + const Matrix & xs, + const Matrix & T, + Tensor3 & Ns) { + for (UInt i = 0; i < natural_coords.cols(); ++i) { + Matrix N_T = Ns(i); + Vector X = natural_coords(i); + Matrix N(interpolation_property::nb_degree_of_freedom, N_T.cols()); + + computeShapes(X, xs, N); + N_T.mul(N.block(0, 0, N_T.rows(), N_T.cols()), T); + } + } + /// compute shape derivatives (input is dxds) for a set of points static inline void computeShapeDerivatives(const Tensor3 & Js, const Tensor3 & DNDSs, const Matrix & R, Tensor3 & Bs) { for (UInt i = 0; i < Js.size(2); ++i) { Matrix J = Js(i); Matrix DNDS = DNDSs(i); Matrix DNDX(DNDS.rows(), DNDS.cols()); auto inv_J = J.inverse(); DNDX.mul(inv_J, DNDS); Matrix B_R = Bs(i); Matrix B(B_R.rows(), B_R.cols()); arrangeInVoigt(DNDX, B); B_R.mul(B, R); } } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with variation of natural coordinates on a given set * of points in natural coordinates */ static inline void computeDNDS(const Matrix & natural_coord, const Matrix & real_coord, Tensor3 & dnds) { for (UInt i = 0; i < natural_coord.cols(); ++i) { Matrix dnds_t = dnds(i); computeDNDS(natural_coord(i), real_coord, dnds_t); } } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with * variation of natural coordinates on a given point in natural * coordinates */ static inline void computeDNDS(const Vector & natural_coord, const Matrix & real_coord, Matrix & dnds); /** * arrange B in Voigt notation from DNDS */ static inline void arrangeInVoigt(const Matrix & dnds, Matrix & B) { // Default implementation assumes dnds is already in Voigt notation B.deepCopy(dnds); } public: static inline constexpr auto getNbNodesPerInterpolationElement() { return interpolation_property::nb_nodes_per_element; } static inline constexpr auto getShapeSize() { return interpolation_property::nb_nodes_per_element * interpolation_property::nb_degree_of_freedom * interpolation_property::nb_degree_of_freedom; } + static inline constexpr auto getShapeIndependantSize() { + return interpolation_property::nb_nodes_per_element * + interpolation_property::nb_degree_of_freedom * + interpolation_property::nb_stress_components; + } static inline constexpr auto getShapeDerivativesSize() { return interpolation_property::nb_nodes_per_element * interpolation_property::nb_degree_of_freedom * interpolation_property::nb_stress_components; } static inline constexpr auto getNaturalSpaceDimension() { return interpolation_property::natural_space_dimension; } static inline constexpr auto getNbDegreeOfFreedom() { return interpolation_property::nb_degree_of_freedom; } static inline constexpr auto getNbStressComponents() { return interpolation_property::nb_stress_components; } }; /// Macro to generate the element class structures for different structural /// element types /* -------------------------------------------------------------------------- */ #define AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY( \ elem_type, geom_type, interp_type, parent_el_type, elem_kind, sp, \ gauss_int_type, min_int_order) \ template <> struct ElementClassProperty { \ static const GeometricalType geometrical_type{geom_type}; \ static const InterpolationType interpolation_type{interp_type}; \ static const ElementType parent_element_type{parent_el_type}; \ static const ElementKind element_kind{elem_kind}; \ static const UInt spatial_dimension{sp}; \ static const GaussIntegrationType gauss_integration_type{gauss_int_type}; \ static const UInt polynomial_degree{min_int_order}; \ } /* -------------------------------------------------------------------------- */ /* ElementClass for structural elements */ /* -------------------------------------------------------------------------- */ template class ElementClass : public GeometricalElement< ElementClassProperty::geometrical_type>, public InterpolationElement< ElementClassProperty::interpolation_type> { protected: using geometrical_element = GeometricalElement::geometrical_type>; using interpolation_element = InterpolationElement< ElementClassProperty::interpolation_type>; using parent_element = ElementClass::parent_element_type>; public: static inline void computeRotationMatrix(Matrix & /*R*/, const Matrix & /*X*/, const Vector & /*extra_normal*/) { AKANTU_TO_IMPLEMENT(); } /// compute jacobian (or integration variable change factor) for a given point static inline void computeJMat(const Vector & natural_coords, const Matrix & Xs, Matrix & J) { Matrix dnds(Xs.rows(), Xs.cols()); parent_element::computeDNDS(natural_coords, dnds); J.mul(dnds, Xs); } static inline void computeJMat(const Matrix & natural_coords, const Matrix & Xs, Tensor3 & Js) { for (UInt i = 0; i < natural_coords.cols(); ++i) { // because non-const l-value reference does not bind to r-value Matrix J = Js(i); computeJMat(Vector(natural_coords(i)), Xs, J); } } static inline void computeJacobian(const Matrix & natural_coords, const Matrix & node_coords, Vector & jacobians) { using itp = typename interpolation_element::interpolation_property; Tensor3 Js(itp::natural_space_dimension, itp::natural_space_dimension, natural_coords.cols()); computeJMat(natural_coords, node_coords, Js); for (UInt i = 0; i < natural_coords.cols(); ++i) { Matrix J = Js(i); jacobians(i) = J.det(); } } static inline void computeRotation(const Matrix & node_coords, Matrix & rotation); public: static AKANTU_GET_MACRO_NOT_CONST(Kind, _ek_structural, ElementKind); static AKANTU_GET_MACRO_NOT_CONST(P1ElementType, _not_defined, ElementType); static AKANTU_GET_MACRO_NOT_CONST(FacetType, _not_defined, ElementType); static constexpr auto getFacetType(__attribute__((unused)) UInt t = 0) { return _not_defined; } static constexpr AKANTU_GET_MACRO_NOT_CONST( SpatialDimension, ElementClassProperty::spatial_dimension, UInt); static constexpr auto getFacetTypes() { return ElementClass<_not_defined>::getFacetTypes(); } }; } // namespace akantu /* -------------------------------------------------------------------------- */ #include "element_class_hermite_inline_impl.hh" /* keep order */ #include "element_class_bernoulli_beam_inline_impl.hh" #include "element_class_kirchhoff_shell_inline_impl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_ELEMENT_CLASS_STRUCTURAL_HH_ */ diff --git a/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.hh b/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.hh index 8145600a2..98b59a223 100644 --- a/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.hh +++ b/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.hh @@ -1,223 +1,233 @@ /** * @file element_class_bernoulli_beam_inline_impl.hh * * @author Fabian Barras * @author Lucas Frerot * * @date creation: Fri Jul 15 2011 * @date last modification: Mon Feb 19 2018 * * @brief Specialization of the element_class class for the type * _bernoulli_beam_2 * * * Copyright (©) 2010-2018 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 . * * * @verbatim --x-----q1----|----q2-----x---> x -1 0 1 @endverbatim * */ /* -------------------------------------------------------------------------- */ #include "aka_static_if.hh" #include "element_class_structural.hh" //#include "aka_element_classes_info.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_CLASS_BERNOULLI_BEAM_INLINE_IMPL_HH_ #define AKANTU_ELEMENT_CLASS_BERNOULLI_BEAM_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY(_itp_bernoulli_beam_2, _itp_lagrange_segment_2, 3, 2, 6); AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY(_itp_bernoulli_beam_3, _itp_lagrange_segment_2, 6, 4, 6); AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY(_bernoulli_beam_2, _gt_segment_2, _itp_bernoulli_beam_2, _segment_2, _ek_structural, 2, _git_segment, 3); AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY(_bernoulli_beam_3, _gt_segment_2, _itp_bernoulli_beam_3, _segment_2, _ek_structural, 3, _git_segment, 3); /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_bernoulli_beam_2, _itk_structural>::computeShapes( const Vector & natural_coords, const Matrix & real_coord, Matrix & N) { Vector L(2); InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeShapes( natural_coords, L); Matrix H(2, 4); InterpolationElement<_itp_hermite_2, _itk_structural>::computeShapes( natural_coords, real_coord, H); // clang-format off // u1 v1 t1 u2 v2 t2 N = {{L(0), 0 , 0 , L(1), 0 , 0 }, // u {0 , H(0, 0), H(0, 1), 0 , H(0, 2), H(0, 3)}, // v {0 , H(1, 0), H(1, 1), 0 , H(1, 2), H(1, 3)}}; // theta // clang-format on } template <> inline void InterpolationElement<_itp_bernoulli_beam_3, _itk_structural>::computeShapes( const Vector & natural_coords, const Matrix & real_coord, Matrix & N) { Vector L(2); InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeShapes( natural_coords, L); Matrix H(2, 4); InterpolationElement<_itp_hermite_2, _itk_structural>::computeShapes( natural_coords, real_coord, H); // clang-format off - // u1 v1 w1 x1 y1 z1 u2 v2 w2 x2 y2 z2 - N = {{L(0), 0 , 0 , 0 , 0 , 0 , L(1), 0 , 0 , 0 , 0 , 0 }, // u - {0 , H(0, 0), 0 , 0 , H(0, 1), 0 , 0 , H(0, 2), 0 , 0 , H(0, 3), 0 }, // v - {0 , 0 , H(0, 0), 0 , 0 , H(0, 1), 0 , 0 , H(0, 2), 0 , 0 , H(0, 3)}, // w - {0 , 0 , 0 , L(0), 0 , 0 , 0 , 0 , 0 , L(1), 0 , 0 }, // thetax - {0 , H(1, 0), 0 , 0 , H(1, 1), 0 , 0 , H(1, 2), 0 , 0 , H(1, 3), 0 }, // thetay - {0 , 0 , H(1, 0), 0 , 0 , H(1, 1), 0 , 0 , H(1, 2), 0 , 0 , H(1, 3)}}; // thetaz + // u1 v1 w1 tx1 ty1 tz1 u2 v2 w2 tx2 ty2 tz2 + N = {{L(0), 0 , 0 , 0 , 0 , 0 , L(1), 0 , 0 , 0 , 0 , 0 }, // u + {0 , H(0, 0), 0 , 0 , 0 , H(0, 1), 0 , H(0, 2), 0 , 0 , 0 , H(0, 3)}, // v + {0 , 0 , H(0, 0), 0 , -H(0, 1), 0 , 0 , 0 , H(0, 2), 0 , -H(0, 3), 0 }, // w + {0 , 0 , 0 , L(0), 0 , 0 , 0 , 0 , 0 , L(1), 0 , 0 }, // thetax + {0 , 0 , H(1, 0), 0 , -H(1, 1), 0 , 0 , 0 , H(1, 2), 0 , -H(1, 3), 0 }, // thetay + {0 , H(1, 0), 0 , 0 , 0 , H(1, 1), 0 , H(1, 2), 0 , 0 , 0 , H(1, 3)}}; // thetaz // clang-format on } +/* -------------------------------------------------------------------------- */ +#if 0 +template <> +inline void +InterpolationElement<_itp_bernoulli_beam_3, _itk_structural>::computeShapesDisplacements( + const Vector & natural_coords, const Matrix & real_coord, + Matrix & N) { +} +#endif + /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_bernoulli_beam_2, _itk_structural>::computeDNDS( const Vector & natural_coords, const Matrix & real_coord, Matrix & dnds) { Matrix L(1, 2); InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeDNDS( natural_coords, L); Matrix H(1, 4); InterpolationElement<_itp_hermite_2, _itk_structural>::computeDNDS( natural_coords, real_coord, H); // Storing the derivatives in dnds dnds.block(L, 0, 0); dnds.block(H, 0, 2); } /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_bernoulli_beam_2, _itk_structural>::arrangeInVoigt( const Matrix & dnds, Matrix & B) { auto L = dnds.block(0, 0, 1, 2); // Lagrange shape derivatives auto H = dnds.block(0, 2, 1, 4); // Hermite shape derivatives // clang-format off - // u1 v1 t1 u2 v2 t2 - B = {{L(0, 0), 0, 0, L(0, 1), 0, 0 }, - {0, H(0, 0), H(0, 1), 0, H(0, 2), H(0, 3)}}; + // u1 v1 t1 u2 v2 t2 + B = {{L(0, 0), 0, 0, L(0, 1), 0, 0 }, + {0, -H(0, 0), -H(0, 1), 0, -H(0, 2), -H(0, 3)}}; // clang-format on } /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_bernoulli_beam_3, _itk_structural>::computeDNDS( const Vector & natural_coords, const Matrix & real_coord, Matrix & dnds) { InterpolationElement<_itp_bernoulli_beam_2, _itk_structural>::computeDNDS( natural_coords, real_coord, dnds); } /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_bernoulli_beam_3, _itk_structural>::arrangeInVoigt( const Matrix & dnds, Matrix & B) { auto L = dnds.block(0, 0, 1, 2); // Lagrange shape derivatives auto H = dnds.block(0, 2, 1, 4); // Hermite shape derivatives // clang-format off - // u1 v1 w1 x1 y1 z1 u2 v2 w2 x2 y2 z2 - B = {{L(0, 0), 0 , 0 , 0 , 0 , 0 , L(0, 1), 0 , 0 , 0 , 0 , 0 }, // eps - {0 , H(0, 0), 0 , 0 , 0 , H(0, 1), 0 , H(0, 2), 0 , 0 , 0 , H(0, 3)}, // chi strong axis - {0 , 0 ,-H(0, 0), 0 , H(0, 1), 0 , 0 , 0 ,-H(0, 2), 0 , H(0, 3), 0 }, // chi weak axis - {0 , 0 , 0 , L(0, 0), 0 , 0 , 0 , 0 , 0 , L(0, 1), 0 , 0 }}; // chi torsion + // u1 v1 w1 x1 y1 z1 u2 v2 w2 x2 y2 z2 + B = {{L(0, 0), 0 , 0 , 0 , 0 , 0 , L(0, 1), 0 , 0 , 0 , 0 , 0 }, // eps + {0 , -H(0, 0), 0 , 0 , 0 , -H(0, 1), 0 , -H(0, 2), 0 , 0 , 0 ,-H(0, 3)}, // chi strong axis + {0 , 0 , -H(0, 0), 0 , H(0, 1) , 0 , 0 , 0 , -H(0, 2) , 0 , H(0, 3) , 0 }, // chi weak axis + {0 , 0 , 0 , L(0, 0), 0 , 0 , 0 , 0 , 0 , L(0, 1), 0 , 0 }}; // chi torsion // clang-format on } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_bernoulli_beam_2>::computeRotationMatrix( Matrix & R, const Matrix & X, const Vector & /*n*/) { Vector x2 = X(1); // X2 Vector x1 = X(0); // X1 auto cs = (x2 - x1); cs.normalize(); auto c = cs(0); auto s = cs(1); // clang-format off /// Definition of the rotation matrix R = {{ c, s, 0.}, {-s, c, 0.}, { 0., 0., 1.}}; // clang-format on } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_bernoulli_beam_3>::computeRotationMatrix( Matrix & R, const Matrix & X, const Vector & n) { Vector x2 = X(1); // X2 Vector x1 = X(0); // X1 auto dim = X.rows(); auto x = (x2 - x1); x.normalize(); auto x_n = x.crossProduct(n); Matrix Pe = {{1., 0., 0.}, {0., -1., 0.}, {0., 0., 1.}}; Matrix Pg(dim, dim); Pg(0) = x; Pg(1) = x_n; Pg(2) = n; Pe *= Pg.inverse(); R.zero(); /// Definition of the rotation matrix for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { R(i + dim, j + dim) = R(i, j) = Pe(i, j); } } } } // namespace akantu #endif /* AKANTU_ELEMENT_CLASS_BERNOULLI_BEAM_INLINE_IMPL_HH_ */ diff --git a/src/fe_engine/element_classes/element_class_hermite_inline_impl.hh b/src/fe_engine/element_classes/element_class_hermite_inline_impl.hh index cbb1c5949..77ad5d0ff 100644 --- a/src/fe_engine/element_classes/element_class_hermite_inline_impl.hh +++ b/src/fe_engine/element_classes/element_class_hermite_inline_impl.hh @@ -1,175 +1,176 @@ /** * @file element_class_hermite_inline_impl.hh * * @author Fabian Barras * @author Lucas Frerot * * @date creation: Fri Nov 10 2017 * @date last modification: Mon Feb 19 2018 * * @brief Specialization of the element_class class for the type * _hermite * * * Copyright (©) 2016-2018 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 . * * * @verbatim --x-----q1----|----q2-----x---> x -1 0 1 @endverbatim * * @f[ * \begin{array}{ll} * M_1(\xi) &= 1/4(\xi^{3}/-3\xi+2)\\ * M_2(\xi) &= -1/4(\xi^{3}-3\xi-2) * \end{array} * * \begin{array}{ll} * L_1(\xi) &= 1/4(\xi^{3}-\xi^{2}-\xi+1)\\ * L_2(\xi) &= 1/4(\xi^{3}+\xi^{2}-\xi-1) * \end{array} * * \begin{array}{ll} * M'_1(\xi) &= 3/4(\xi^{2}-1)\\ * M'_2(\xi) &= -3/4(\xi^{2}-1) * \end{array} * * \begin{array}{ll} * L'_1(\xi) &= 1/4(3\xi^{2}-2\xi-1)\\ * L'_2(\xi) &= 1/4(3\xi^{2}+2\xi-1) * \end{array} *@f] * * *@f[ * \begin{array}{ll} * N'_1(\xi) &= -1/2\\ * N'_2(\xi) &= 1/2 * \end{array}] * * \begin{array}{ll} * -M''_1(\xi) &= -3\xi/2\\ * -M''_2(\xi) &= 3\xi/2\\ * \end{array} * * \begin{array}{ll} * -L''_1(\xi) &= -1/2a(3\xi/a-1)\\ * -L''_2(\xi) &= -1/2a(3\xi/a+1) * \end{array} *@f] * */ /* -------------------------------------------------------------------------- */ #include "aka_static_if.hh" #include "element_class_structural.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_CLASS_HERMITE_INLINE_IMPL_HH_ #define AKANTU_ELEMENT_CLASS_HERMITE_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY(_itp_hermite_2, _itp_lagrange_segment_2, 2, 1, 4); /* -------------------------------------------------------------------------- */ namespace { namespace details { inline Real computeLength(const Matrix & real_coord) { Vector x1 = real_coord(0); Vector x2 = real_coord(1); return x1.distance(x2); } inline void computeShapes(const Vector & natural_coords, Real a, Matrix & N) { /// natural coordinate Real xi = natural_coords(0); - + auto xi2 = xi * xi; + auto xi3 = xi * xi * xi; // Cubic Hermite splines interpolating displacement - auto M1 = 1. / 4. * Math::pow<2>(xi - 1) * (xi + 2); - auto M2 = 1. / 4. * Math::pow<2>(xi + 1) * (2 - xi); - auto L1 = a / 4. * Math::pow<2>(xi - 1) * (xi + 1); - auto L2 = a / 4. * Math::pow<2>(xi + 1) * (xi - 1); + auto M1 = 1. / 4. * (2. - 3. * xi + xi3); + auto M2 = 1. / 4. * (2. + 3. * xi - xi3); + auto L1 = a / 4. * (1 - xi - xi2 + xi3); + auto L2 = a / 4. * (-1 - xi + xi2 + xi3);; #if 1 // Version where we also interpolate the rotations // Derivatives (with respect to x) of previous functions interpolating // rotations - auto M1_ = 3. / (4. * a) * (xi * xi - 1); - auto M2_ = 3. / (4. * a) * (1 - xi * xi); - auto L1_ = 1 / 4. * (3 * xi * xi - 2 * xi - 1); - auto L2_ = 1 / 4. * (3 * xi * xi + 2 * xi - 1); + auto M1_ = 3. / (4. * a) * (xi2 - 1); + auto M2_ = 3. / (4. * a) * (1 - xi2); + auto L1_ = 1 / 4. * (3 * xi2 - 2 * xi - 1); + auto L2_ = 1 / 4. * (3 * xi2 + 2 * xi - 1); // clang-format off // v1 t1 v2 t2 N = {{M1 , L1 , M2 , L2}, // displacement interpolation {M1_, L1_, M2_, L2_}}; // rotation interpolation // clang-format on #else // Version where we only interpolate displacements // clang-format off // v1 t1 v2 t2 N = {{M1, L1, M2, L2}}; // clang-format on #endif } /* ---------------------------------------------------------------------- */ inline void computeDNDS(const Vector & natural_coords, Real a, Matrix & B) { // natural coordinate Real xi = natural_coords(0); // Derivatives with respect to xi for rotations auto M1 = 3. / 2. * xi; auto M2 = 3. / 2. * (-xi); - auto L1 = a / 2. * (3 * xi - 1); - auto L2 = a / 2. * (3 * xi + 1); + auto L1 = 1. * a / 2. * (3 * xi - 1); + auto L2 = 1. * a / 2. * (3 * xi + 1); // v1 t1 v2 t2 B = {{M1, L1, M2, L2}}; // computing curvature : {chi} = [B]{d} B /= a; // to account for first order deriv w/r to x } } // namespace details } // namespace /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_hermite_2, _itk_structural>::computeShapes( const Vector & natural_coords, const Matrix & real_coord, Matrix & N) { auto L = details::computeLength(real_coord); details::computeShapes(natural_coords, L / 2, N); } /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_hermite_2, _itk_structural>::computeDNDS( const Vector & natural_coords, const Matrix & real_coord, Matrix & B) { auto L = details::computeLength(real_coord); details::computeDNDS(natural_coords, L / 2, B); } } // namespace akantu #endif /* AKANTU_ELEMENT_CLASS_HERMITE_INLINE_IMPL_HH_ */ diff --git a/src/fe_engine/fe_engine_template_tmpl_field.hh b/src/fe_engine/fe_engine_template_tmpl_field.hh index c217c5892..bd1b4d5fb 100644 --- a/src/fe_engine/fe_engine_template_tmpl_field.hh +++ b/src/fe_engine/fe_engine_template_tmpl_field.hh @@ -1,455 +1,503 @@ /** * @file fe_engine_template_tmpl_field.hh * * @author Nicolas Richart * * @date creation: Wed Aug 09 2017 * @date last modification: Thu Dec 07 2017 * * @brief implementation of the assemble field s functions * * * Copyright (©) 2016-2018 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 "fe_engine_template.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_FE_ENGINE_TEMPLATE_TMPL_FIELD_HH_ #define AKANTU_FE_ENGINE_TEMPLATE_TMPL_FIELD_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ /* Matrix lumping functions */ /* -------------------------------------------------------------------------- */ namespace fe_engine { namespace details { namespace { template void fillField(const Functor & field_funct, Array & field, UInt nb_element, UInt nb_integration_points, ElementType type, GhostType ghost_type) { UInt nb_degree_of_freedom = field.getNbComponent(); field.resize(nb_integration_points * nb_element); auto field_it = field.begin_reinterpret( nb_degree_of_freedom, nb_integration_points, nb_element); Element el{type, 0, ghost_type}; for (; el.element < nb_element; ++el.element, ++field_it) { field_funct(*field_it, el); } } } // namespace } // namespace details } // namespace fe_engine /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template struct AssembleLumpedTemplateHelper { template