diff --git a/examples/python/custom-material-dynamics/material.dat b/examples/python/custom-material-dynamics/material.dat deleted file mode 100644 index cbb2e7b68..000000000 --- a/examples/python/custom-material-dynamics/material.dat +++ /dev/null @@ -1,4 +0,0 @@ -material local_elastic [ - name = steel - rho = 1 # density -] diff --git a/examples/python/custom-material-dynamics/newmark.py b/examples/python/custom-material-dynamics/newmark.py deleted file mode 100644 index 351eab73f..000000000 --- a/examples/python/custom-material-dynamics/newmark.py +++ /dev/null @@ -1,180 +0,0 @@ -#!/usr/bin/python -################################################################ -import os -import subprocess -import numpy as np -import akantu -################################################################ - -class FixedValue: - - def __init__(self, value, axis): - self.value = value - if axis == 'x': - axis = 0 - if axis == 'y': - axis = 1 - self.axis = axis - - def operator(self, node, flags, disp, coord): - # sets the displacement to the desired value in the desired axis - disp[self.axis] = self.value - # sets the blocked dofs vector to true in the desired axis - flags[self.axis] = True - -################################################################ - - -class LocalElastic: - - def __init__(self): - - ## young modulus - self.E = 1 - ## Poisson coefficient - self.nu = 0.3 - ## density - self.rho = 1 - ## First Lame coefficient - self._lambda = self.nu * self.E / ((1 + self.nu) * (1 - 2*self.nu)) - ## Second Lame coefficient (shear modulus) - self.mu = self.E / (2 * (1 + self.nu)) - - ## declares all the internals - @staticmethod - def registerInternals(): - return [] - - ## declares all the parameters that could be parsed - @staticmethod - def registerParam(): - return [] - - ## declares all the parameters that are needed - def getPushWaveSpeed(self): - return np.sqrt((self._lambda + 2*self.mu)/self.rho) - - ## compute small deformation tensor - @staticmethod - def computeEpsilon(grad_u): - return 0.5 * (grad_u + np.einsum('aij->aji', grad_u)) - - ## constitutive law for a given quadrature point - def computeStress(self, grad_u,sigma, dummy_internals): - lbda = 1. - mu = 1. - nquads = grad_u.shape[0] - grad_u = grad_u.reshape((nquads, 2, 2)) - epsilon = self.computeEpsilon(grad_u) - sigma = sigma.reshape((nquads, 2, 2)) - trace = np.trace(grad_u, axis1=1, axis2=2) - sigma[:, :, :] = (np.einsum('a,ij->aij', trace, lbda*np.eye(2)) - + mu * epsilon) - - -################################################################ -def main(): - - spatial_dimension = 2 - Lbar = 10. - akantu.initialize('material.dat') - - mesh_file = 'bar.msh' - max_steps = 250 - time_step = 1e-3 - - #if mesh was not created the calls gmsh to generate it - if not os.path.isfile(mesh_file): - ret = subprocess.call('gmsh -2 bar.geo bar.msh', shell=True) - if ret != 0: - raise Exception('execution of GMSH failed: do you have it installed ?') - - - ################################################################ - ## Initialization - ################################################################ - mesh = akantu.Mesh(spatial_dimension) - mesh.read(mesh_file) - - mesh.createGroupsFromStringMeshData("physical_names") - model = akantu.SolidMechanicsModel(mesh) - - model.initFull(akantu.SolidMechanicsModelOptions(akantu._explicit_lumped_mass, True)) - mat = LocalElastic() - model.registerNewPythonMaterial(mat, "local_elastic") - model.initMaterials() - - - model.setBaseName("waves") - model.addDumpFieldVector("displacement") - model.addDumpFieldVector("acceleration") - model.addDumpFieldVector("velocity") - model.addDumpFieldVector("residual") - model.addDumpFieldVector("force") - model.addDumpField("strain") - model.addDumpField("stress") - model.addDumpField("blocked_dofs") - - ################################################################ - ## boundary conditions - ################################################################ - - model.applyDirichletBC(FixedValue(0, 'x'), "XBlocked") - model.applyDirichletBC(FixedValue(0, 'y'), "YBlocked") - - ################################################################ - ## initial conditions - ################################################################ - - displacement = model.getDisplacement() - nb_nodes = mesh.getNbNodes() - position = mesh.getNodes() - - pulse_width = 1 - A = 0.01 - for i in range(0, nb_nodes): - # Sinus * Gaussian - x = position[i, 0] - 5. - L = pulse_width - k = 0.1 * 2 * np.pi * 3 / L - displacement[i, 0] = A * np.sin(k * x) * np.exp(-(k * x) * (k * x) / (L * L)) - - ################################################################ - ## timestep value computation - ################################################################ - time_factor = 0.8 - stable_time_step = model.getStableTimeStep() * time_factor - - print "Stable Time Step = {0}".format(stable_time_step) - print "Required Time Step = {0}".format(time_step) - - time_step = stable_time_step * time_factor - - model.setTimeStep(time_step) - - ################################################################ - ## loop for evolution of motion dynamics - ################################################################ - model.updateResidual() - - epot = model.getEnergy('potential') - ekin = model.getEnergy('kinetic') - - print "step,step * time_step,epot,ekin,epot + ekin" - for step in range(0, max_steps+1): - - if step % 10 == 0: - model.dump() - ## output energy calculation to screen - print "{0},{1},{2},{3},{4}".format(step, step * time_step, - epot, ekin, - (epot + ekin)) - - model.solveStep() - - akantu.finalize() - return - -################################################################ -if __name__ == "__main__": - main() diff --git a/examples/python/custom-material-dynamics/CMakeLists.txt b/examples/python/custom-material/CMakeLists.txt similarity index 100% rename from examples/python/custom-material-dynamics/CMakeLists.txt rename to examples/python/custom-material/CMakeLists.txt diff --git a/examples/python/custom-material-dynamics/bar.geo b/examples/python/custom-material/bar.geo similarity index 95% copy from examples/python/custom-material-dynamics/bar.geo copy to examples/python/custom-material/bar.geo index 32f173f35..dc70c3402 100644 --- a/examples/python/custom-material-dynamics/bar.geo +++ b/examples/python/custom-material/bar.geo @@ -1,45 +1,44 @@ // Mesh size h = 0.05; h1 = h; h2 = h; // Dimensions of the bar Lx = 10; Ly = 1; // ------------------------------------------ // Geometry // ------------------------------------------ Point(101) = { 0.0, -Ly/2, 0.0, h1}; Point(102) = { Lx, -Ly/2, 0.0, h2}; Point(103) = { Lx, 0., 0.0, h2}; Point(104) = { Lx, Ly/2., 0.0, h2}; Point(105) = { 0.0, Ly/2., 0.0, h1}; Point(106) = { 0.0, 0., 0.0, h1}; Line(101) = {101,102}; Line(102) = {102,103}; Line(103) = {103,104}; Line(104) = {104,105}; Line(105) = {105,106}; Line(106) = {106,101}; Line(107) = {106,103}; Line Loop(108) = {101, 102, -107, 106}; Plane Surface(109) = {108}; Line Loop(110) = {103, 104, 105, 107}; Plane Surface(111) = {110}; Physical Surface(112) = {109, 111}; Transfinite Surface "*"; Recombine Surface "*"; Physical Surface(113) = {111, 109}; Physical Line("XBlocked") = {103, 102}; -Physical Line("ImposedVelocity") = {105, 106}; Physical Line("YBlocked") = {104, 101}; diff --git a/examples/python/custom-material/bi-material.py b/examples/python/custom-material/bi-material.py new file mode 100644 index 000000000..017acc202 --- /dev/null +++ b/examples/python/custom-material/bi-material.py @@ -0,0 +1,192 @@ +from __future__ import print_function +# ------------------------------------------------------------- # +import akantu +import subprocess +import numpy as np +import time +# ------------------------------------------------------------- # + + +class LocalElastic: + + # declares all the internals + def initMaterial(self, internals, params): + self.E = params['E'] + self.nu = params['nu'] + self.rho = params['rho'] + # print(self.__dict__) + # First Lame coefficient + self.lame_lambda = self.nu * self.E / ( + (1. + self.nu) * (1. - 2. * self.nu)) + # Second Lame coefficient (shear modulus) + self.lame_mu = self.E / (2. * (1. + self.nu)) + + all_factor = internals['factor'] + all_quad_coords = internals['quad_coordinates'] + + for elem_type in all_factor.keys(): + factor = all_factor[elem_type] + quad_coords = all_quad_coords[elem_type] + + factor[:] = 1. + factor[quad_coords[:, 1] < 0.5] = .5 + + # declares all the internals + @staticmethod + def registerInternals(): + return ['potential', 'factor'] + + # declares all the internals + @staticmethod + def registerInternalSizes(): + return [1, 1] + + # declares all the parameters that could be parsed + @staticmethod + def registerParam(): + return ['E', 'nu'] + + # declares all the parameters that are needed + def getPushWaveSpeed(self, params): + return np.sqrt((self.lame_lambda + 2 * self.lame_mu) / self.rho) + + # compute small deformation tensor + @staticmethod + def computeEpsilon(grad_u): + return 0.5 * (grad_u + np.einsum('aij->aji', grad_u)) + + # constitutive law + def computeStress(self, grad_u, sigma, internals, params): + n_quads = grad_u.shape[0] + grad_u = grad_u.reshape((n_quads, 2, 2)) + factor = internals['factor'].reshape(n_quads) + epsilon = self.computeEpsilon(grad_u) + sigma = sigma.reshape((n_quads, 2, 2)) + trace = np.einsum('aii->a', grad_u) + + sigma[:, :, :] = ( + np.einsum('a,ij->aij', trace, + self.lame_lambda * np.eye(2)) + + 2. * self.lame_mu * epsilon) + + # print(sigma.reshape((n_quads, 4))) + # print(grad_u.reshape((n_quads, 4))) + sigma[:, :, :] = np.einsum('aij, a->aij', sigma, factor) + + # constitutive law tangent modulii + def computeTangentModuli(self, grad_u, tangent, internals, params): + n_quads = tangent.shape[0] + tangent = tangent.reshape(n_quads, 3, 3) + factor = internals['factor'].reshape(n_quads) + + Miiii = self.lame_lambda + 2 * self.lame_mu + Miijj = self.lame_lambda + Mijij = self.lame_mu + + tangent[:, 0, 0] = Miiii + tangent[:, 1, 1] = Miiii + tangent[:, 0, 1] = Miijj + tangent[:, 1, 0] = Miijj + tangent[:, 2, 2] = Mijij + tangent[:, :, :] = np.einsum('aij, a->aij', tangent, factor) + + # computes the energy density + def getEnergyDensity(self, energy_type, energy_density, + grad_u, stress, internals, params): + + nquads = stress.shape[0] + stress = stress.reshape(nquads, 2, 2) + grad_u = grad_u.reshape((nquads, 2, 2)) + + if energy_type != 'potential': + raise RuntimeError('not known energy') + + epsilon = self.computeEpsilon(grad_u) + + energy_density[:, 0] = ( + 0.5 * np.einsum('aij,aij->a', stress, epsilon)) + + +# applies manually the boundary conditions +def applyBC(model): + + nbNodes = model.getMesh().getNbNodes() + position = model.getMesh().getNodes() + displacement = model.getDisplacement() + blocked_dofs = model.getBlockedDOFs() + + width = 1. + height = 1. + epsilon = 1e-8 + for node in range(0, nbNodes): + if((np.abs(position[node, 0]) < epsilon) or # left side + (np.abs(position[node, 0] - width) < epsilon)): # right side + blocked_dofs[node, 0] = True + displacement[node, 0] = 0 * position[node, 0] + 0. + + if(np.abs(position[node, 1]) < epsilon): # lower side + blocked_dofs[node, 1] = True + displacement[node, 1] = - 1. + if(np.abs(position[node, 1] - height) < epsilon): # upper side + blocked_dofs[node, 1] = True + displacement[node, 1] = 1. + + +################################################################ +# main +################################################################ + +# main paraeters +spatial_dimension = 2 +mesh_file = 'square.msh' + +# call gmsh to generate the mesh +ret = subprocess.call(['gmsh', '-2', 'square.geo', '-optimize', 'square.msh']) +if ret != 0: + raise Exception( + 'execution of GMSH failed: do you have it installed ?') + +time.sleep(1) + +# read mesh +mesh = akantu.Mesh(spatial_dimension) +mesh.read(mesh_file) + +# create the custom material +mat = LocalElastic() +akantu.registerNewPythonMaterial(mat, "local_elastic") + +# init akantu +akantu.initialize('material.dat') + +# init the SolidMechanicsModel +model = akantu.SolidMechanicsModel(mesh) +model.initFull(_analysis_method=akantu._static) + +# configure the solver +solver = model.getNonLinearSolver() +solver.set("max_iterations", 2) +solver.set("threshold", 1e-3) +solver.set("convergence_type", akantu._scc_solution) + +# prepare the dumper +model.setBaseName("bimaterial") +model.addDumpFieldVector("displacement") +model.addDumpFieldVector("internal_force") +model.addDumpFieldVector("external_force") +model.addDumpField("strain") +model.addDumpField("stress") +model.addDumpField("factor") +model.addDumpField("blocked_dofs") + +# Boundary conditions +applyBC(model) + +# solve the problem +model.solveStep() + +# dump paraview files +model.dump() + +# shutdown akantu +akantu.finalize() diff --git a/examples/python/custom-material/custom-material.py b/examples/python/custom-material/custom-material.py new file mode 100644 index 000000000..199d51dbe --- /dev/null +++ b/examples/python/custom-material/custom-material.py @@ -0,0 +1,210 @@ +#!/usr/bin/env python3 + +from __future__ import print_function +################################################################ +import os +import subprocess +import numpy as np +import akantu +################################################################ + + +class FixedValue: + + def __init__(self, value, axis): + self.value = value + self.axis = axis + + def operator(self, node, flags, disp, coord): + # sets the displacement to the desired value in the desired axis + disp[self.axis] = self.value + # sets the blocked dofs vector to true in the desired axis + flags[self.axis] = True + +################################################################ + + +class LocalElastic: + + # declares all the internals + def initMaterial(self, internals, params): + self.E = params['E'] + self.nu = params['nu'] + self.rho = params['rho'] + # First Lame coefficient + self.lame_lambda = self.nu * self.E / ( + (1 + self.nu) * (1 - 2 * self.nu)) + # Second Lame coefficient (shear modulus) + self.lame_mu = self.E / (2 * (1 + self.nu)) + + # declares all the internals + @staticmethod + def registerInternals(): + return ['potential'] + + # declares all the internals + @staticmethod + def registerInternalSizes(): + return [1] + + # declares all the parameters that could be parsed + @staticmethod + def registerParam(): + return ['E', 'nu'] + + # declares all the parameters that are needed + def getPushWaveSpeed(self, params): + return np.sqrt((self.lame_lambda + 2 * self.lame_mu) / self.rho) + + # compute small deformation tensor + @staticmethod + def computeEpsilon(grad_u): + return 0.5 * (grad_u + np.einsum('aij->aji', grad_u)) + + # constitutive law + def computeStress(self, grad_u, sigma, internals, params): + nquads = grad_u.shape[0] + grad_u = grad_u.reshape((nquads, 2, 2)) + epsilon = self.computeEpsilon(grad_u) + sigma = sigma.reshape((nquads, 2, 2)) + trace = np.trace(grad_u, axis1=1, axis2=2) + sigma[:, :, :] = ( + np.einsum('a,ij->aij', trace, + self.lame_lambda * np.eye(2)) + + 2.*self.lame_mu * epsilon) + + # constitutive law tangent modulii + def computeTangentModuli(self, grad_u, tangent, internals, params): + n_quads = tangent.shape[0] + tangent = tangent.reshape(n_quads, 3, 3) + + Miiii = self.lame_lambda + 2 * self.lame_mu + Miijj = self.lame_lambda + Mijij = self.lame_mu + + tangent[:, 0, 0] = Miiii + tangent[:, 1, 1] = Miiii + tangent[:, 0, 1] = Miijj + tangent[:, 1, 0] = Miijj + tangent[:, 2, 2] = Mijij + + # computes the energy density + def getEnergyDensity(self, energy_type, energy_density, + grad_u, stress, internals, params): + + nquads = stress.shape[0] + stress = stress.reshape(nquads, 2, 2) + grad_u = grad_u.reshape((nquads, 2, 2)) + + if energy_type != 'potential': + raise RuntimeError('not known energy') + + epsilon = self.computeEpsilon(grad_u) + + energy_density[:, 0] = ( + 0.5 * np.einsum('aij,aij->a', stress, epsilon)) + + +################################################################ +# main +################################################################ + +spatial_dimension = 2 +akantu.initialize('material.dat') + +mesh_file = 'bar.msh' +max_steps = 250 +time_step = 1e-3 + +# if mesh was not created the calls gmsh to generate it +if not os.path.isfile(mesh_file): + ret = subprocess.call('gmsh -2 bar.geo bar.msh', shell=True) + if ret != 0: + raise Exception( + 'execution of GMSH failed: do you have it installed ?') + +################################################################ +# Initialization +################################################################ +mesh = akantu.Mesh(spatial_dimension) +mesh.read(mesh_file) + +mat = LocalElastic() +akantu.registerNewPythonMaterial(mat, "local_elastic") + +model = akantu.SolidMechanicsModel(mesh) + +model.initFull(_analysis_method=akantu._explicit_lumped_mass) +# model.initFull(_analysis_method=akantu._implicit_dynamic) + +model.setBaseName("waves") +model.addDumpFieldVector("displacement") +model.addDumpFieldVector("acceleration") +model.addDumpFieldVector("velocity") +model.addDumpFieldVector("internal_force") +model.addDumpFieldVector("external_force") +model.addDumpField("strain") +model.addDumpField("stress") +model.addDumpField("blocked_dofs") + +################################################################ +# boundary conditions +################################################################ + +model.applyDirichletBC(FixedValue(0, akantu._x), "XBlocked") +model.applyDirichletBC(FixedValue(0, akantu._y), "YBlocked") + +################################################################ +# initial conditions +################################################################ + +displacement = model.getDisplacement() +nb_nodes = mesh.getNbNodes() +position = mesh.getNodes() + +pulse_width = 1 +A = 0.01 +for i in range(0, nb_nodes): + # Sinus * Gaussian + x = position[i, 0] - 5. + L = pulse_width + k = 0.1 * 2 * np.pi * 3 / L + displacement[i, 0] = A * \ + np.sin(k * x) * np.exp(-(k * x) * (k * x) / (L * L)) + +################################################################ +# timestep value computation +################################################################ +time_factor = 0.8 +stable_time_step = model.getStableTimeStep() * time_factor + +print("Stable Time Step = {0}".format(stable_time_step)) +print("Required Time Step = {0}".format(time_step)) + +time_step = stable_time_step * time_factor + +model.setTimeStep(time_step) + +################################################################ +# loop for evolution of motion dynamics +################################################################ +model.assembleInternalForces() + +print("step,step * time_step,epot,ekin,epot + ekin") +for step in range(0, max_steps + 1): + + model.solveStep() + + if step % 10 == 0: + model.dump() + + epot = model.getEnergy('potential') + ekin = model.getEnergy('kinetic') + + # output energy calculation to screen + print("{0},{1},{2},{3},{4}".format(step, step * time_step, + epot, ekin, + (epot + ekin))) + +akantu.finalize() + diff --git a/examples/python/custom-material/material.dat b/examples/python/custom-material/material.dat new file mode 100644 index 000000000..1498fc3c9 --- /dev/null +++ b/examples/python/custom-material/material.dat @@ -0,0 +1,6 @@ +material local_elastic [ + name = fictive + rho = 1. # density + E = 1. # young modulus + nu = .3 # poisson ratio +] diff --git a/examples/python/custom-material/square.geo b/examples/python/custom-material/square.geo new file mode 100644 index 000000000..120adaf42 --- /dev/null +++ b/examples/python/custom-material/square.geo @@ -0,0 +1,33 @@ +// Mesh size +h = .05; + +h1 = h; +h2 = h; + +// Dimensions of the bar +Lx = 1.; +Ly = 1.; + +// ------------------------------------------ +// Geometry +// ------------------------------------------ + +Point(101) = { 0.0, 0.0, 0.0, h1}; +Point(102) = { Lx, 0.0, 0.0, h2}; + +Point(103) = { Lx, Ly, 0.0, h2}; +Point(104) = { 0., Ly, 0.0, h2}; +//+ +Line(1) = {101, 102}; +//+ +Line(2) = {102, 103}; +//+ +Line(3) = {103, 104}; +//+ +Line(4) = {104, 101}; +//+ +Line Loop(1) = {4, 1, 2, 3}; +//+ +Plane Surface(1) = {1}; +//+ +Physical Surface("bulk") = {1}; diff --git a/examples/python/custom-material-dynamics/bar.geo b/examples/python/dynamics/bar.geo similarity index 100% rename from examples/python/custom-material-dynamics/bar.geo rename to examples/python/dynamics/bar.geo diff --git a/examples/python/dynamics/dynamics.py b/examples/python/dynamics/dynamics.py new file mode 100644 index 000000000..23d34898e --- /dev/null +++ b/examples/python/dynamics/dynamics.py @@ -0,0 +1,130 @@ +#!/usr/bin/env python3 + +from __future__ import print_function +################################################################ +import os +import subprocess +import numpy as np +import akantu +################################################################ + + +class FixedValue: + + def __init__(self, value, axis): + self.value = value + self.axis = axis + + def operator(self, node, flags, disp, coord): + # sets the displacement to the desired value in the desired axis + disp[self.axis] = self.value + # sets the blocked dofs vector to true in the desired axis + flags[self.axis] = True + +################################################################ + + +def main(): + + spatial_dimension = 2 + + akantu.initialize('material.dat') + + mesh_file = 'bar.msh' + max_steps = 250 + time_step = 1e-3 + + # if mesh was not created the calls gmsh to generate it + if not os.path.isfile(mesh_file): + ret = subprocess.call('gmsh -2 bar.geo bar.msh', shell=True) + if ret != 0: + raise Exception( + 'execution of GMSH failed: do you have it installed ?') + + ################################################################ + # Initialization + ################################################################ + mesh = akantu.Mesh(spatial_dimension) + mesh.read(mesh_file) + + model = akantu.SolidMechanicsModel(mesh) + + model.initFull(_analysis_method=akantu._explicit_lumped_mass) + # model.initFull(_analysis_method=akantu._implicit_dynamic) + + model.setBaseName("waves") + model.addDumpFieldVector("displacement") + model.addDumpFieldVector("acceleration") + model.addDumpFieldVector("velocity") + model.addDumpFieldVector("internal_force") + model.addDumpFieldVector("external_force") + model.addDumpField("strain") + model.addDumpField("stress") + model.addDumpField("blocked_dofs") + + ################################################################ + # boundary conditions + ################################################################ + + model.applyDirichletBC(FixedValue(0, akantu._x), "XBlocked") + model.applyDirichletBC(FixedValue(0, akantu._y), "YBlocked") + + ################################################################ + # initial conditions + ################################################################ + + displacement = model.getDisplacement() + nb_nodes = mesh.getNbNodes() + position = mesh.getNodes() + + pulse_width = 1 + A = 0.01 + for i in range(0, nb_nodes): + # Sinus * Gaussian + x = position[i, 0] - 5. + L = pulse_width + k = 0.1 * 2 * np.pi * 3 / L + displacement[i, 0] = A * \ + np.sin(k * x) * np.exp(-(k * x) * (k * x) / (L * L)) + + ################################################################ + # timestep value computation + ################################################################ + time_factor = 0.8 + stable_time_step = model.getStableTimeStep() * time_factor + + print("Stable Time Step = {0}".format(stable_time_step)) + print("Required Time Step = {0}".format(time_step)) + + time_step = stable_time_step * time_factor + + model.setTimeStep(time_step) + + ################################################################ + # loop for evolution of motion dynamics + ################################################################ + model.assembleInternalForces() + + print("step,step * time_step,epot,ekin,epot + ekin") + for step in range(0, max_steps + 1): + + model.solveStep() + + if step % 10 == 0: + model.dump() + + epot = model.getEnergy('potential') + ekin = model.getEnergy('kinetic') + + # output energy calculation to screen + print("{0},{1},{2},{3},{4}".format(step, step * time_step, + epot, ekin, + (epot + ekin))) + + akantu.finalize() + return + + +################################################################ +if __name__ == "__main__": + main() diff --git a/examples/python/dynamics/material.dat b/examples/python/dynamics/material.dat new file mode 100644 index 000000000..42b9b8ae7 --- /dev/null +++ b/examples/python/dynamics/material.dat @@ -0,0 +1,6 @@ +material elastic [ + name = fictive + rho = 1. # density + E = 1. # young modulus + nu = .3 # poisson ratio +] diff --git a/examples/python/plate-hole/plate.py b/examples/python/plate-hole/plate.py index d02f6e181..bc821c0ea 100644 --- a/examples/python/plate-hole/plate.py +++ b/examples/python/plate-hole/plate.py @@ -1,107 +1,113 @@ +#!/usr/bin/env python3 + +from __future__ import print_function + import akantu import numpy as np ################################################################ -## Dirichlet Boudary condition functor: fix the displacement +# Dirichlet Boudary condition functor: fix the displacement ################################################################ + + class FixedValue: - def __init__(self,value,axis): + def __init__(self, value, axis): self.value = value - if axis == 'x': axis = 0 - if axis == 'y': axis = 1 self.axis = axis - - def operator(self,node,flags,disp,coord): + + def operator(self, node, flags, disp, coord): # sets the displacement to the desired value in the desired axis disp[self.axis] = self.value # sets the blocked dofs vector to true in the desired axis flags[self.axis] = True ################################################################ -## Neumann Boudary condition functor: impose a traction +# Neumann Boudary condition functor: impose a traction ################################################################ + + class FromTraction: - def __init__(self,traction): + def __init__(self, traction): self.traction = traction - def operator(self,quad_point,force,coord,normals): + def operator(self, quad_point, force, coord, normals): # sets the force to the desired value in the desired axis force[:] = self.traction ################################################################ -def solve(material_file,mesh_file,traction): + +def solve(material_file, mesh_file, traction): akantu.initialize(material_file) spatial_dimension = 2 - + ################################################################ - ##Initialization + # Initialization ################################################################ mesh = akantu.Mesh(spatial_dimension) mesh.read(mesh_file) - akantu.MeshUtils.purifyMesh(mesh) - mesh.createGroupsFromStringMeshData("physical_names") model = akantu.SolidMechanicsModel(mesh) model.initFull(akantu.SolidMechanicsModelOptions(akantu._static)) - model.assembleStiffnessMatrix() - model.updateResidual() - + model.setBaseName("plate") model.addDumpFieldVector("displacement") - model.addDumpFieldVector("force") - model.addDumpField("boundary") + model.addDumpFieldVector("external_force") model.addDumpField("strain") model.addDumpField("stress") model.addDumpField("blocked_dofs") - + ################################################################ - ##Boundary conditions + # Boundary conditions ################################################################ - displacement = model.getDisplacement() - blocked_dofs = model.getBlockedDOFs() - - model.applyDirichletBC(FixedValue(0.0, 'x'), "XBlocked") - model.applyDirichletBC(FixedValue(0.0, 'y'), "YBlocked") + + model.applyDirichletBC(FixedValue(0.0, akantu._x), "XBlocked") + model.applyDirichletBC(FixedValue(0.0, akantu._y), "YBlocked") trac = np.zeros(spatial_dimension) trac[1] = traction - print "Solve for traction ",traction + print("Solve for traction ", traction) model.getForce()[:] = 0 model.applyNeumannBC(FromTraction(trac), "Traction") - model.solveStaticDisplacement(1e-10,2); + solver = model.getNonLinearSolver() + solver.set("max_iterations", 2) + solver.set("threshold", 1e-10) + solver.set("convergence_type", akantu._scc_residual) + + model.solveStep() model.dump() akantu.finalize() ################################################################ -# main +# main ################################################################ + def main(): import os mesh_file = 'plate.msh' - #if mesh was not created the calls gmsh to generate it + # if mesh was not created the calls gmsh to generate it if not os.path.isfile(mesh_file): import subprocess - ret = subprocess.call('gmsh -2 plate.geo {0}'.format(mesh_file),shell=True) + ret = subprocess.call( + 'gmsh -2 plate.geo {0}'.format(mesh_file), shell=True) if not ret == 0: - raise Exception('execution of GMSH failed: do you have it installed ?') + raise Exception( + 'execution of GMSH failed: do you have it installed ?') - material_file = 'material.dat' - spatial_dimension = 2 traction = 1. - solve(material_file,mesh_file,traction) + solve(material_file, mesh_file, traction) + ################################################################ if __name__ == "__main__": - main() - + main() diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 335cce7eb..f42b5fb03 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -1,271 +1,271 @@ #=============================================================================== # @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 . # #=============================================================================== #=============================================================================== # Configuration #=============================================================================== package_get_all_definitions(AKANTU_DEFS) list(REMOVE_ITEM AKANTU_DEFS AKANTU_CORE_CXX11) #message(${AKANTU_DEFS}) set(AKA_DEFS "") foreach (def ${AKANTU_DEFS}) list(APPEND AKA_DEFS "-D${def}") endforeach() -set(AKANTU_SWIG_FLAGS -w309,325,401,317,509,503,383,384 ${AKA_DEFS}) +set(AKANTU_SWIG_FLAGS -w309,325,401,317,509,503,383,384,476,362 ${AKA_DEFS}) set(AKANTU_SWIG_OUTDIR ${CMAKE_CURRENT_SOURCE_DIR}) set(AKANTU_SWIG_MODULES swig/akantu.i) #=============================================================================== # Swig wrapper #=============================================================================== set(SWIG_REQURIED_VERISON 3.0) find_package(SWIG ${SWIG_REQURIED_VERISON}) find_package(PythonInterp ${AKANTU_PREFERRED_PYTHON_VERSION} REQUIRED) mark_as_advanced(SWIG_EXECUTABLE) if(NOT PYTHON_VERSION_MAJOR LESS 3) list(APPEND AKANTU_SWIG_FLAGS -py3) endif() package_get_all_include_directories( AKANTU_LIBRARY_INCLUDE_DIRS ) set(_swig_include_dirs ${CMAKE_CURRENT_SOURCE_DIR}/swig ${AKANTU_LIBRARY_INCLUDE_DIRS} ${PROJECT_BINARY_DIR}/src ${AKANTU_EXTERNAL_INCLUDE_DIR} ) include(CMakeParseArguments) function(swig_generate_dependencies _module _depedencies _depedencies_out) set(_dependencies_script "${CMAKE_CURRENT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/_swig_generate_dependencies.cmake") file(WRITE ${_dependencies_script} " set(_include_directories ${_include_directories}) list(APPEND _include_directories \"./\") set(_dep) set(_files_to_process \${_module}) while(_files_to_process) list(GET _files_to_process 0 _file) list(REMOVE_AT _files_to_process 0) file(STRINGS \${_file} _file_content REGEX \"^%include *\\\"(.*)\\\"\") set(_includes) foreach(_line \${_file_content}) string(REGEX REPLACE \"^%include *\\\"(.*)\\\"\" \"\\\\1\" _inc \${_line}) if(_inc) list(APPEND _includes \${_inc}) endif() endforeach() foreach(_include \${_includes}) unset(_found) foreach(_inc_dir \${_include_directories}) if(EXISTS \${_inc_dir}/\${_include}) set(_found \${_inc_dir}/\${_include}) break() endif() endforeach() if(_found) list(APPEND _files_to_process \${_found}) list(APPEND _dep \${_found}) endif() endforeach() endwhile() get_filename_component(_module_we \"\${_module}\" NAME_WE) set(_dependencies_file \${CMAKE_CURRENT_BINARY_DIR}\${CMAKE_FILES_DIRECTORY}/_swig_\${_module_we}_depends.cmake) file(WRITE \"\${_dependencies_file}\" \"set(_swig_\${_module_we}_depends\") foreach(_d \${_dep}) file(APPEND \"\${_dependencies_file}\" \" \${_d}\") endforeach() file(APPEND \"\${_dependencies_file}\" \" )\") ") get_filename_component(_module_absolute "${_module}" ABSOLUTE) get_filename_component(_module_we "${_module}" NAME_WE) set(_dependencies_file ${CMAKE_CURRENT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/_swig_${_module_we}_depends.cmake) if(EXISTS ${_dependencies_file}) include(${_dependencies_file}) else() execute_process(COMMAND ${CMAKE_COMMAND} -D_module=${_module_absolute} -P ${_dependencies_script} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) include(${_dependencies_file}) endif() set(${_depedencies_out} ${_swig_${_module_we}_depends} PARENT_SCOPE) add_custom_command(OUTPUT ${_dependencies_file} COMMAND ${CMAKE_COMMAND} -D_module=${_module_absolute} -P ${_dependencies_script} COMMENT "Scanning dependencies for swig module ${_module_we}" WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} MAIN_DEPENDENCY ${_module_absolute} DEPENDS ${_swig_${_module_we}_depends} ) set(${_depedencies} ${_dependencies_file} PARENT_SCOPE) endfunction() function(swig_generate_wrappers project _wrappers_cpp _wrappers_py) cmake_parse_arguments(_swig_opt "" "OUTPUT_DIR;DEPENDENCIES" "EXTRA_FLAGS;INCLUDE_DIRECTORIES" ${ARGN}) if(_swig_opt_OUTPUT_DIR) set(_output_dir ${_swig_opt_OUTPUT_DIR}) else() set(_output_dir ${CMAKE_CURRENT_BINARY_DIR}) endif() set(_swig_wrappers) get_directory_property(_include_directories INCLUDE_DIRECTORIES) list(APPEND _include_directories ${_swig_opt_INCLUDE_DIRECTORIES}) if(_include_directories) string(REPLACE ";" ";-I" _swig_include_directories "${_include_directories}") endif() foreach(_module ${_swig_opt_UNPARSED_ARGUMENTS}) swig_generate_dependencies(${_module} _module_dependencies _depends_out) if(_swig_opt_DEPENDENCIES) set(${_swig_opt_DEPENDENCIES} ${_depends_out} PARENT_SCOPE) endif() get_filename_component(_module_absolute "${_module}" ABSOLUTE) get_filename_component(_module_path "${_module_absolute}" PATH) get_filename_component(_module_name "${_module}" NAME) get_filename_component(_module_we "${_module}" NAME_WE) set(_wrapper "${_output_dir}/${_module_we}_wrapper.cpp") set(_extra_wrapper "${_output_dir}/${_module_we}.py") set(_extra_wrapper_bin "${CMAKE_CURRENT_BINARY_DIR}/${_module_we}.py") if(SWIG_FOUND) set_source_files_properties("${_wrapper}" PROPERTIES GENERATED 1) set_source_files_properties("${_extra_wrapper}" PROPERTIES GENERATED 1) set(_dependencies_file ${CMAKE_CURRENT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/_swig_${_module_we}_depends.cmake) set(_ouput "${_wrapper}" "${_extra_wrapper}") add_custom_command( OUTPUT ${_ouput} COMMAND "${SWIG_EXECUTABLE}" ARGS -python -c++ ${_swig_opt_EXTRA_FLAGS} -outdir ${_output_dir} -I${_swig_include_directories} -I${_module_path} -o "${_wrapper}" "${_module_absolute}" COMMAND ${CMAKE_COMMAND} -E copy_if_different ${_extra_wrapper} ${_extra_wrapper_bin} # MAIN_DEPENDENCY "${_module_absolute}" DEPENDS ${_module_dependencies} COMMENT "Generating swig wrapper ${_module} -> ${_wrapper}" ) list(APPEND _swig_wrappers ${_wrapper}) list(APPEND _swig_wrappers_py "${_extra_wrapper_bin}") else() if(NOT EXISTS ${_wrapper} OR NOT EXISTS "${_extra_wrapper}") message(FATAL_ERROR "The file ${_wrapper} and/or ${_extra_wrapper} does " "not exists and they cannot be generated. Install swig ${SWIG_REQURIED_VERISON} " " in order to generate them. Or get them from a different machine, " "in order to be able to compile the python interface") else() list(APPEND _swig_wrappers "${_wrapper}") list(APPEND _swig_wrappers_py "${_extra_wrapper_bin}") endif() endif() endforeach() add_custom_target(${project}_generate_swig_wrappers DEPENDS ${_swig_wrappers}) set(${_wrappers_cpp} ${_swig_wrappers} PARENT_SCOPE) set(${_wrappers_py} ${_swig_wrappers_py} PARENT_SCOPE) endfunction() swig_generate_wrappers(akantu AKANTU_SWIG_WRAPPERS_CPP AKANTU_WRAPPERS_PYTHON ${AKANTU_SWIG_MODULES} EXTRA_FLAGS ${AKANTU_SWIG_FLAGS} DEPENDENCIES _deps INCLUDE_DIRECTORIES ${_swig_include_dirs}) if(AKANTU_SWIG_WRAPPERS_CPP) set(_ext_files "${AKANTU_SWIG_WRAPPERS_CPP}") set(_inc_dirs "${_swig_include_dirs}") set(_lib_dirs "${Akantu_BINARY_DIR}/src") string(TOUPPER "${CMAKE_BUILD_TYPE}" _config) set(_akantu_lib_name "akantu${CMAKE_${_config}_POSTFIX}") get_property(_compile_flags TARGET akantu PROPERTY COMPILE_FLAGS) set(_compile_flags "${_compile_flags} ${CMAKE_CXX_FLAGS} ${CMAKE_CXX_FLAGS_${_config}}") separate_arguments(_compile_flags) get_property(_cxx_standard TARGET akantu PROPERTY CXX_STANDARD) set(_flags ${_compile_flags} -std=c++${_cxx_standard} -Wno-maybe-uninitialized -Wno-unused-parameter) set(_quiet) if(NOT CMAKE_VERBOSE_MAKEFILE) set(_quiet --quiet) endif() configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in ${CMAKE_CURRENT_BINARY_DIR}/setup.py @ONLY) add_custom_target(_akantu ALL COMMAND ${PYTHON_EXECUTABLE} ./setup.py ${_quiet} --no-user-cfg build_ext --inplace WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} DEPENDS ${AKANTU_SWIG_WRAPPERS_CPP} akantu COMMENT "Building akantu's python interface" ) set_directory_properties(PROPERTIES ADDITIONAL_MAKE_CLEAN_FILES ${CMAKE_CURRENT_BINARY_DIR}/_akantu${CMAKE_SHARED_MODULE_SUFFIX}) install(CODE "execute_process( COMMAND ${PYTHON_EXECUTABLE} ./setup.py ${_quiet} install --prefix=${AKANTU_PYTHON_INSTALL_PREFIX} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})" ) package_declare_extra_files_to_package(python_interface ${_deps} ${PROJECT_SOURCE_DIR}/python/${AKANTU_SWIG_MODULES}) endif() diff --git a/python/swig/aka_common.i b/python/swig/aka_common.i index 51fc2b460..fe3a0f149 100644 --- a/python/swig/aka_common.i +++ b/python/swig/aka_common.i @@ -1,124 +1,142 @@ /** * @file aka_common.i * * @author Guillaume Anciaux * @author Fabian Barras * @author Nicolas Richart * * @date creation: Fri Dec 12 2014 * @date last modification: Wed Jan 13 2016 * * @brief wrapper to aka_common.hh * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ %{ //#include "aka_config.hh" #include "aka_common.hh" #include "aka_csr.hh" #include "element.hh" +#include "python_functor.hh" %} +%include "stl.i" + namespace akantu { %ignore getStaticParser; %ignore getUserParser; %ignore ghost_types; %ignore initialize(int & argc, char ** & argv); %ignore initialize(const std::string & input_file, int & argc, char ** & argv); extern const Array empty_filter; } %typemap(in) (int argc, char *argv[]) { int i = 0; if (!PyList_Check($input)) { PyErr_SetString(PyExc_ValueError, "Expecting a list"); return NULL; } $1 = PyList_Size($input); $2 = new char *[$1+1]; for (i = 0; i < $1; i++) { PyObject *s = PyList_GetItem($input,i); if (!PyString_Check(s)) { free($2); PyErr_SetString(PyExc_ValueError, "List items must be strings"); return NULL; } $2[i] = PyString_AsString(s); } $2[i] = 0; } %typemap(freearg) (int argc, char *argv[]) { %#if defined(__INTEL_COMPILER) //#pragma warning ( disable : 383 ) %#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu %#elif (defined(__GNUC__) || defined(__GNUG__)) %# if __cplusplus > 199711L %# pragma GCC diagnostic ignored "-Wmaybe-uninitialized" %# endif %#endif delete [] $2; %#if defined(__INTEL_COMPILER) //#pragma warning ( disable : 383 ) %#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu %#elif (defined(__GNUC__) || defined(__GNUG__)) %# if __cplusplus > 199711L %# pragma GCC diagnostic pop %# endif %#endif } %inline %{ namespace akantu { #if defined(AKANTU_USE_MPI) const int MPI=1; #else const int MPI=0; #endif void _initializeWithoutArgv(const std::string & input_file) { int nb_args = 0; char ** null = nullptr; initialize(input_file, nb_args, null); } + +#define AKANTU_PP_STR_TO_TYPE2(s, data, elem) \ + ({ BOOST_PP_STRINGIZE(elem) , elem}) + + PyObject * getElementTypes(){ + + std::map element_types{ + BOOST_PP_SEQ_FOR_EACH_I( + AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(AKANTU_ek_regular_ELEMENT_TYPE), + BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE2, akantu, AKANTU_ek_regular_ELEMENT_TYPE))}; + + return akantu::PythonFunctor::convertToPython(element_types); + } } + %} %pythoncode %{ import sys as _aka_sys def initialize(input_file="", argv=_aka_sys.argv): if _aka_sys.modules[__name__].MPI == 1: try: from mpi4py import MPI except ImportError: pass _initializeWithoutArgv(input_file) %} %include "aka_config.hh" %include "aka_common.hh" %include "aka_element_classes_info.hh" %include "element.hh" +%include "aka_error.hh" diff --git a/python/swig/akantu.i b/python/swig/akantu.i index 59ce7182e..aecdc6d80 100644 --- a/python/swig/akantu.i +++ b/python/swig/akantu.i @@ -1,82 +1,80 @@ /** * @file akantu.i * * @author Guillaume Anciaux * @author Fabian Barras * @author Nicolas Richart * * @date creation: Fri Dec 12 2014 * @date last modification: Mon Nov 23 2015 * * @brief Main swig file for akantu' python interface * * @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 . * */ %module akantu %exception { try { $action } catch (akantu::debug::Exception e) { PyErr_SetString(PyExc_IndexError,e.what()); return NULL; } } -%include "stl.i" - #define __attribute__(x) %ignore akantu::operator <<; %include "aka_common.i" %include "aka_csr.i" %include "aka_array.i" %define print_self(MY_CLASS) %extend akantu::MY_CLASS { std::string __str__() { std::stringstream sstr; sstr << *($self); return sstr.str(); } } %enddef %include "mesh.i" %include "mesh_utils.i" %include "model.i" %include "communicator.i" %include "solid_mechanics_model.i" #if defined(AKANTU_COHESIVE_ELEMENT) %include "solid_mechanics_model_cohesive.i" #endif #if defined(AKANTU_HEAT_TRANSFER) %include "heat_transfer_model.i" #endif #if defined(AKANTU_STRUCTURAL_MECHANICS) %include "load_functions.i" %include "structural_mechanics_model.i" #endif diff --git a/python/swig/heat_transfer_model.i b/python/swig/heat_transfer_model.i index e4398ba93..c6c61a1a8 100644 --- a/python/swig/heat_transfer_model.i +++ b/python/swig/heat_transfer_model.i @@ -1,52 +1,79 @@ /** * @file heat_transfer_model.i * * @author Guillaume Anciaux * * @date creation: Wed Jul 15 2015 * * @brief heat transfer model wrapper * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ %{ #include "heat_transfer_model.hh" #include "data_accessor.hh" + #include "python_functor.hh" %} namespace akantu { %ignore HeatTransferModel::initFEEngineBoundary; %ignore HeatTransferModel::initParallel; %ignore HeatTransferModel::initArrays; %ignore HeatTransferModel::initMaterials; %ignore HeatTransferModel::initModel; %ignore HeatTransferModel::initPBC; %ignore HeatTransferModel::initSolver; %ignore HeatTransferModel::getNbDataToPack; %ignore HeatTransferModel::getNbData; %ignore HeatTransferModel::packData; %ignore HeatTransferModel::unpackData; } %include "heat_transfer_model.hh" + + +%extend akantu::HeatTransferModel { + + Real getParamReal(const ID & param){ + Real res = $self->get(param); + return res; + } + UInt getParamUInt(const ID & param){ + UInt res = $self->get(param); + return res; + } + int getParamInt(const ID & param){ + int res = $self->get(param); + return res; + } + + PyObject * getParamMatrix(const ID & param){ + Matrix res = $self->get(param); + // I know it is ugly !!!!! sorry nico: missing time to do something clean + Matrix * res2 = new Matrix(res); + PyObject * pobj = akantu::PythonFunctor::convertToPython(*res2); + return pobj; + } +} + diff --git a/python/swig/material.i b/python/swig/material.i index c1d3b152c..9d8d4c948 100644 --- a/python/swig/material.i +++ b/python/swig/material.i @@ -1,50 +1,75 @@ /** * @file material.i * * @author Nicolas Richart * * @brief material wrapper * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ %{ #include "solid_mechanics_model.hh" #include "material_python.hh" -%} + #include "parameter_registry.hh" + #include "parsable.hh" + #include "parser.hh" + %} namespace akantu { %ignore Material::onNodesAdded; %ignore Material::onNodesRemoved; %ignore Material::onElementsAdded; %ignore Material::onElementsRemoved; %ignore Material::onElementsChanged; %ignore Material::getParam; } %include "material.hh" %extend akantu::Material { Array & getArrayReal(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) { return $self->getArray(id, type, ghost_type); } -} + + Real getParamReal(const ID & param){ + Real res = $self->get(param); + return res; + } + UInt getParamUInt(const ID & param){ + UInt res = $self->get(param); + return res; + } + int getParamInt(const ID & param){ + int res = $self->get(param); + return res; + } + + PyObject * getParamMatrix(const ID & param){ + Matrix res = $self->get(param); + // I know it is ugly !!!!! sorry nico: missing time to do something clean + Matrix * res2 = new Matrix(res); + PyObject * pobj = akantu::PythonFunctor::convertToPython(*res2); + return pobj; + } + } + diff --git a/python/swig/mesh.i b/python/swig/mesh.i index 5711ced7b..08016e7f9 100644 --- a/python/swig/mesh.i +++ b/python/swig/mesh.i @@ -1,160 +1,166 @@ /** * @file mesh.i * * @author Guillaume Anciaux * @author Fabian Barras * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Dec 12 2014 * @date last modification: Wed Jan 13 2016 * * @brief mesh wrapper * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ %{ #include "mesh.hh" #include "node_group.hh" #include "solid_mechanics_model.hh" -// #include "dumpable_inline_impl.hh" - +#include "python_functor.hh" +#include "mesh_utils.hh" + using akantu::IntegrationPoint; using akantu::Vector; using akantu::ElementTypeMapArray; using akantu::MatrixProxy; using akantu::Matrix; using akantu::UInt; using akantu::Real; using akantu::Array; using akantu::SolidMechanicsModel; %} namespace akantu { %ignore NewNodesEvent; %ignore RemovedNodesEvent; %ignore NewElementsEvent; %ignore RemovedElementsEvent; %ignore MeshEventHandler; %ignore MeshEvent< UInt >; %ignore MeshEvent< Element >; %ignore Mesh::extractNodalCoordinatesFromPBCElement; %ignore Mesh::getGroupDumer; %ignore Mesh::getFacetLocalConnectivity; %ignore Mesh::getAllFacetTypes; %ignore Mesh::getCommunicator; %ignore GroupManager::getElementGroups; %ignore Dumpable::addDumpFieldExternalReal; } print_self(Mesh) // Swig considers enums to be ints, and it creates a conflict with two versions of getNbElement() %rename(getNbElementByDimension) akantu::Mesh::getNbElement(const UInt spatial_dimension = _all_dimensions, const GhostType& ghost_type = _not_ghost, const ElementKind& kind = _ek_not_defined) const; %extend akantu::Mesh { + PyObject * getElementGroups(){ + return akantu::PythonFunctor::convertToPython($self->getElementGroups()); + } + void resizeMesh(UInt nb_nodes, UInt nb_element, const ElementType & type) { Array & nodes = const_cast &>($self->getNodes()); nodes.resize(nb_nodes); $self->addConnectivityType(type); Array & connectivity = const_cast &>($self->getConnectivity(type)); connectivity.resize(nb_element); } } %extend akantu::GroupManager { void createGroupsFromStringMeshData(const std::string & dataset_name) { if (dataset_name == "physical_names"){ AKANTU_EXCEPTION("Deprecated behavior: no need to call 'createGroupsFromStringMeshData' for physical names"); } $self->createGroupsFromMeshData(dataset_name); } void createGroupsFromUIntMeshData(const std::string & dataset_name) { $self->createGroupsFromMeshData(dataset_name); } } %extend akantu::NodeGroup { akantu::Array & getGroupedNodes(akantu::Array & surface_array, Mesh & mesh) { akantu::Array group_node = $self->getNodes(); akantu::Array & full_array = mesh.getNodes(); surface_array.resize(group_node.size()); for (UInt i = 0; i < group_node.size(); ++i) { for (UInt cmp = 0; cmp < full_array.getNbComponent(); ++cmp) { surface_array(i,cmp) = full_array(group_node(i),cmp); } } akantu::Array & res(surface_array); return res; } akantu::Array & getGroupedArray(akantu::Array & surface_array, akantu::SolidMechanicsModel & model, int type) { akantu::Array * full_array; switch (type) { case 0 : full_array = new akantu::Array(model.getDisplacement()); break; case 1 : full_array = new akantu::Array(model.getVelocity()); break; case 2 : full_array = new akantu::Array(model.getForce()); break; } akantu::Array group_node = $self->getNodes(); surface_array.resize(group_node.size()); for (UInt i = 0; i < group_node.size(); ++i) { for (UInt cmp = 0; cmp < full_array->getNbComponent(); ++cmp) { surface_array(i,cmp) = (*full_array)(group_node(i),cmp); } } akantu::Array & res(surface_array); return res; } } %include "group_manager.hh" -%include "element_group.hh" %include "node_group.hh" %include "dumper_iohelper.hh" %include "dumpable_iohelper.hh" +%include "element_group.hh" %include "mesh.hh" +%include "mesh_utils.hh" namespace akantu{ %extend Dumpable { void addDumpFieldExternalReal(const std::string & field_id, const Array & field){ $self->addDumpFieldExternal(field_id,field); } } } diff --git a/python/swig/model.i b/python/swig/model.i index 93ab9b994..365a13d2b 100644 --- a/python/swig/model.i +++ b/python/swig/model.i @@ -1,109 +1,131 @@ /** * @file model.i * * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Dec 12 2014 * @date last modification: Wed Nov 11 2015 * * @brief model wrapper * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ %{ #include "boundary_condition_python_functor.hh" #include "model_solver.hh" #include "non_linear_solver.hh" %} namespace akantu { %ignore Model::createSynchronizerRegistry; %ignore Model::getSynchronizerRegistry; %ignore Model::createParallelSynch; %ignore Model::getDOFSynchronizer; %ignore Model::registerFEEngineObject; %ignore Model::unregisterFEEngineObject; %ignore Model::getFEEngineBoundary; // %ignore Model::getFEEngine; %ignore Model::getFEEngineClass; %ignore Model::getFEEngineClassBoundary; %ignore Model::setParser; %ignore Model::updateDataForNonLocalCriterion; %ignore IntegrationPoint::operator=; %ignore FEEngine::getNbIntegrationPoints; %ignore FEEngine::getShapes; %ignore FEEngine::getShapesDerivatives; %ignore FEEngine::getIntegrationPoints; %ignore FEEngine::getIGFEMElementTypes; %ignore FEEngine::interpolateOnIntegrationPoints(const Array &,ElementTypeMapArray &,const ElementTypeMapArray *) const; %ignore FEEngine::interpolateOnIntegrationPoints(const Array &,ElementTypeMapArray &) const; %ignore FEEngine::interpolateOnIntegrationPoints(const Array &,Array &,UInt,const ElementType&,const GhostType &,const Array< UInt > &) const; %ignore FEEngine::interpolateOnIntegrationPoints(const Array &,Array &,UInt,const ElementType&,const GhostType &) const; %ignore FEEngine::onNodesAdded; %ignore FEEngine::onNodesRemoved; %ignore FEEngine::onElementsAdded; %ignore FEEngine::onElementsChanged; %ignore FEEngine::onElementsRemoved; %ignore FEEngine::elementTypes; } %include "sparse_matrix.i" %include "fe_engine.hh" %rename(FreeBoundaryDirichlet) akantu::BC::Dirichlet::FreeBoundary; %rename(FreeBoundaryNeumann) akantu::BC::Neumann::FreeBoundary; %rename(PythonBoundary) akantu::BC::Dirichlet::PythonFunctor; %include "boundary_condition_functor.hh" %include "boundary_condition.hh" %include "boundary_condition_python_functor.hh" %include "communication_buffer.hh" %include "data_accessor.hh" //%include "synchronizer.hh" //%include "synchronizer_registry.hh" %include "model.hh" %include "non_linear_solver.hh" %extend akantu::Model { + void initFullImpl( + const akantu::ModelOptions & options = akantu::ModelOptions()){ + $self->initFull(options); + }; + + %insert("python") %{ + def initFull(self, *args, **kwargs): + if len(args) == 0: + import importlib as __aka_importlib + _module = __aka_importlib.import_module(self.__module__) + _cls = getattr(_module, "{0}Options".format(self.__class__.__name__)) + options = _cls() + if len(kwargs) > 0: + for key, val in kwargs.items(): + if key[0] == '_': + key = key[1:] + setattr(options, key, val) + else: + options = args[0] + + self.initFullImpl(options) + %} + void solveStep(){ $self->solveStep(); } akantu::NonLinearSolver & getNonLinearSolver(){ return $self->getNonLinearSolver(); } } %extend akantu::NonLinearSolver { - void set(const std::string & id, akantu::Real val){ if (id == "max_iterations") $self->set(id, int(val)); else if (id == "convergence_type") $self->set(id, akantu::SolveConvergenceCriteria(UInt(val))); else $self->set(id, val); } } diff --git a/python/swig/solid_mechanics_model.i b/python/swig/solid_mechanics_model.i index 784c481e9..b4c8f2c31 100644 --- a/python/swig/solid_mechanics_model.i +++ b/python/swig/solid_mechanics_model.i @@ -1,212 +1,204 @@ /** * @file solid_mechanics_model.i * * @author Guillaume Anciaux * @author Fabian Barras * @author Nicolas Richart * * @date creation: Fri Dec 12 2014 * @date last modification: Wed Jan 06 2016 * * @brief solid mechanics model wrapper * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ %{ #include "model_options.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "boundary_condition.hh" #include "boundary_condition_functor.hh" #include "boundary_condition_python_functor.hh" #include "material_selector.hh" #include "material_python.hh" %} namespace akantu { %ignore SolidMechanicsModel::initFEEngineBoundary; %ignore SolidMechanicsModel::initParallel; %ignore SolidMechanicsModel::initArrays; %ignore SolidMechanicsModel::initModel; %ignore SolidMechanicsModel::initPBC; %ignore SolidMechanicsModel::initExplicit; %ignore SolidMechanicsModel::isExplicit; %ignore SolidMechanicsModel::updateCurrentPosition; %ignore SolidMechanicsModel::updateAcceleration; %ignore SolidMechanicsModel::updateIncrement; %ignore SolidMechanicsModel::updatePreviousDisplacement; %ignore SolidMechanicsModel::saveStressAndStrainBeforeDamage; %ignore SolidMechanicsModel::updateEnergiesAfterDamage; %ignore SolidMechanicsModel::solveLumped; %ignore SolidMechanicsModel::explicitPred; %ignore SolidMechanicsModel::explicitCorr; %ignore SolidMechanicsModel::initSolver; %ignore SolidMechanicsModel::initImplicit; %ignore SolidMechanicsModel::initialAcceleration; %ignore SolidMechanicsModel::testConvergence; %ignore SolidMechanicsModel::testConvergenceIncrement; %ignore SolidMechanicsModel::testConvergenceResidual; %ignore SolidMechanicsModel::initVelocityDampingMatrix; %ignore SolidMechanicsModel::getNbData; %ignore SolidMechanicsModel::packData; %ignore SolidMechanicsModel::unpackData; %ignore SolidMechanicsModel::setMaterialSelector; %ignore SolidMechanicsModel::getSolver; %ignore SolidMechanicsModel::getSynchronizer; %ignore Dumpable::registerExternalDumper; %ignore Material::onNodesAdded; %ignore Material::onNodesRemoved; %ignore Material::onElementsAdded; %ignore Material::onElementsRemoved; %ignore Material::onElementsChanged; %template(SolidMechanicsBoundaryCondition) BoundaryCondition; } %include "dumpable.hh" print_self(SolidMechanicsModel) %include "material.i" %include "model_options.hh" %include "solid_mechanics_model.hh" %inline %{ namespace akantu{ void registerNewPythonMaterial(PyObject * obj, const akantu::ID & mat_type) { MaterialFactory::getInstance().registerAllocator( mat_type, [obj](UInt, const ID &, SolidMechanicsModel & model, const ID & id) -> std::unique_ptr { return std::make_unique(model, obj, id); } ); } } %} %extend akantu::SolidMechanicsModel { - - void initFull( - const akantu::SolidMechanicsModelOptions & options = akantu::SolidMechanicsModelOptions()){ - $self->initFull(options); - }; - - /* ------------------------------------------------------------------------ */ void setPhysicalNamesMaterialSelector(){ AKANTU_EXCEPTION("API change: This method needs to be fixed."); /* akantu::MeshDataMaterialSelector * selector = new */ /* akantu::MeshDataMaterialSelector("physical_names", *self); */ /* self->setMaterialSelector(*selector); */ } /* ------------------------------------------------------------------------ */ void getResidual() { AKANTU_EXCEPTION("Deprecated function. You should replace:\n" "model.getResidual()\n" "with\n" "model.getInternalForce()\n"); } void registerNewPythonMaterial(PyObject * obj, const akantu::ID & mat_type) { AKANTU_EXCEPTION("Deprecated function. You should replace:\n" "model.registerNewPythonMaterial(obj, mat_id)\n" "with\n" "akantu.registerNewPythonMaterial(obj, mat_id)\n\n" "This MUST be done before instanciating the model."); /* std::pair */ /* sub_sect = akantu::getStaticParser().getSubSections(akantu::_st_material); */ /* akantu::Parser::const_section_iterator it = sub_sect.first; */ /* for (; it != sub_sect.second; ++it) { */ /* if (it->getName() == mat_type) { */ /* AKANTU_DEBUG_ASSERT($self->materials_names_to_id.find(mat_type) == */ /* $self->materials_names_to_id.end(), */ /* "A material with this name '" */ /* << mat_type << "' has already been registered. " */ /* << "Please use unique names for materials"); */ /* UInt mat_count = $self->materials.size(); */ /* $self->materials_names_to_id[mat_type] = mat_count; */ /* std::stringstream sstr_mat; */ /* sstr_mat << $self->getID() << ":" << mat_count << ":" << mat_type; */ /* akantu::ID mat_id = sstr_mat.str(); */ /* akantu::Material * material = new akantu::MaterialPython(*$self, obj, mat_id); */ /* $self->materials.push_back(material); */ /* material->parseSection(*it); */ /* } */ /* } */ } /* ------------------------------------------------------------------------ */ void updateResidual() { AKANTU_EXCEPTION("Deprecated function. You should replace:\n" "model.updateResidual()\n" "with\n" "model.assembleInternalForces()\n\n" "beware that the total nodal force is your responsability to compute" " since it is now handled by the solver."); } /* ------------------------------------------------------------------------ */ void applyDirichletBC(PyObject * func_obj, const std::string & group_name) { akantu::BC::PythonFunctorDirichlet functor(func_obj); $self->applyBC(functor, group_name); } /* ------------------------------------------------------------------------ */ void applyNeumannBC(PyObject * func_obj, const std::string & group_name) { akantu::BC::PythonFunctorNeumann functor(func_obj); $self->applyBC(functor, group_name); } /* ------------------------------------------------------------------------ */ void applyUniformPressure(Real pressure, const std::string surface_name) { UInt spatial_dimension = $self->getSpatialDimension(); akantu::Matrix surface_stress(spatial_dimension, spatial_dimension, 0.0); for (UInt i = 0; i < spatial_dimension; ++i) { surface_stress(i, i) = -pressure; } $self->applyBC(akantu::BC::Neumann::FromStress(surface_stress), surface_name); } /* ------------------------------------------------------------------------ */ void blockDOF(const std::string surface_name, SpacialDirection direction) { $self->applyBC(akantu::BC::Dirichlet::FixedValue(0.0, direction), surface_name); } } - diff --git a/src/io/parser/parsable.cc b/src/io/parser/parsable.cc index ce512bb24..89b41a458 100644 --- a/src/io/parser/parsable.cc +++ b/src/io/parser/parsable.cc @@ -1,108 +1,109 @@ /** * @file parsable.cc * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Thu Nov 19 2015 * * @brief Parsable implementation * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "parsable.hh" #include "aka_random_generator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Parsable::Parsable(const ParserType & section_type, const ID & id) : section_type(section_type), pid(id) { this->consisder_sub = false; } /* -------------------------------------------------------------------------- */ Parsable::~Parsable() = default; /* -------------------------------------------------------------------------- */ void Parsable::registerSubSection(const ParserType & type, const std::string & name, Parsable & sub_section) { SubSectionKey key(type, name); sub_sections[key] = &sub_section; this->registerSubRegistry(name, sub_section); } /* -------------------------------------------------------------------------- */ void Parsable::parseParam(const ParserParameter & in_param) { auto it = params.find(in_param.getName()); if (it == params.end()) { if (Parser::isPermissive()) { AKANTU_DEBUG_WARNING("No parameter named " << in_param.getName() << " registered in " << pid << "."); return; } else AKANTU_EXCEPTION("No parameter named " << in_param.getName() << " registered in " << pid << "."); } Parameter & param = *(it->second); param.setAuto(in_param); } /* -------------------------------------------------------------------------- */ void Parsable::parseSection(const ParserSection & section) { if (section_type != section.getType()) AKANTU_EXCEPTION("The object " << pid << " is meant to parse section of type " << section_type << ", so it cannot parse section of type " << section.getType()); auto params = section.getParameters(); auto it = params.first; for (; it != params.second; ++it) { parseParam(*it); } auto sit = section.getSubSections().first; for (; sit != section.getSubSections().second; ++sit) { parseSubSection(*sit); } } /* -------------------------------------------------------------------------- */ void Parsable::parseSubSection(const ParserSection & section) { SubSectionKey key(section.getType(), section.getName()); auto it = sub_sections.find(key); if (it != sub_sections.end()) { it->second->parseSection(section); } else if (!Parser::isPermissive()) { AKANTU_EXCEPTION("No parsable defined for sub sections of type <" << key.first << "," << key.second << "> in " << pid); - } else + } else { AKANTU_DEBUG_WARNING("No parsable defined for sub sections of type <" << key.first << "," << key.second << "> in " << pid); + } } } // akantu diff --git a/src/mesh/group_manager.hh b/src/mesh/group_manager.hh index 9a11829c6..1a203623a 100644 --- a/src/mesh/group_manager.hh +++ b/src/mesh/group_manager.hh @@ -1,308 +1,317 @@ /** * @file group_manager.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Mon Nov 16 2015 * * @brief Stores information relevent to the notion of element and nodes *groups. * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GROUP_MANAGER_HH__ #define __AKANTU_GROUP_MANAGER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class ElementGroup; class NodeGroup; class Mesh; class Element; class ElementSynchronizer; template class CommunicationBufferTemplated; namespace dumper { class Field; } } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class GroupManager { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: + +#ifdef SWIGPYTHON +public: using ElementGroups = std::map; using NodeGroups = std::map; +private: +#else + using ElementGroups = std::map; + using NodeGroups = std::map; +#endif + public: using GroupManagerTypeSet = std::set; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GroupManager(const Mesh & mesh, const ID & id = "group_manager", const MemoryID & memory_id = 0); virtual ~GroupManager(); /* ------------------------------------------------------------------------ */ /* Groups iterators */ /* ------------------------------------------------------------------------ */ public: using node_group_iterator = NodeGroups::iterator; using element_group_iterator = ElementGroups::iterator; using const_node_group_iterator = NodeGroups::const_iterator; using const_element_group_iterator = ElementGroups::const_iterator; #ifndef SWIG #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(group_type, function, \ param_in, param_out) \ inline BOOST_PP_CAT(BOOST_PP_CAT(const_, group_type), _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) const { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ }; \ \ inline BOOST_PP_CAT(group_type, _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ } #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(group_type, function) \ AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION( \ group_type, function, BOOST_PP_EMPTY(), BOOST_PP_EMPTY()) AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(element_group, find, const std::string & name, name); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(node_group, find, const std::string & name, name); #endif /* ------------------------------------------------------------------------ */ /* Clustering filter */ /* ------------------------------------------------------------------------ */ public: class ClusteringFilter { public: virtual bool operator()(const Element &) const { return true; } }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// create an empty node group NodeGroup & createNodeGroup(const std::string & group_name, bool replace_group = false); /// create a node group from another node group but filtered template NodeGroup & createFilteredNodeGroup(const std::string & group_name, const NodeGroup & node_group, T & filter); /// destroy a node group void destroyNodeGroup(const std::string & group_name); /// create an element group and the associated node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension = _all_dimensions, bool replace_group = false); /// create an element group from another element group but filtered template ElementGroup & createFilteredElementGroup(const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter); /// destroy an element group and the associated node group void destroyElementGroup(const std::string & group_name, bool destroy_node_group = false); /// destroy all element groups and the associated node groups void destroyAllElementGroups(bool destroy_node_groups = false); /// create a element group using an existing node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group); /// create groups based on values stored in a given mesh data template void createGroupsFromMeshData(const std::string & dataset_name); /// create boundaries group by a clustering algorithm \todo extend to parallel UInt createBoundaryGroupFromGeometry(); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, Mesh & mesh_facets, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); private: /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, const std::string & cluster_name_prefix, const ClusteringFilter & filter, Mesh & mesh_facets); public: /// Create an ElementGroup based on a NodeGroup void createElementGroupFromNodeGroup(const std::string & name, const std::string & node_group, UInt dimension = _all_dimensions); virtual void printself(std::ostream & stream, int indent = 0) const; /// this function insure that the group names are present on all processors /// /!\ it is a SMP call void synchronizeGroupNames(); /// register an elemental field to the given group name (overloading for /// ElementalPartionField) #ifndef SWIG template class dump_type> dumper::Field * createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// ElementalField) template class ret_type, template class, bool> class dump_type> dumper::Field * createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// MaterialInternalField) template class dump_type> dumper::Field * createElementalField(const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); template class ftype> dumper::Field * createNodalField(const ftype * field, const std::string & group_name, UInt padding_size = 0); template class ftype> dumper::Field * createStridedNodalField(const ftype * field, const std::string & group_name, UInt size, UInt stride, UInt padding_size); protected: /// fill a buffer with all the group names void fillBufferWithGroupNames( CommunicationBufferTemplated & comm_buffer) const; /// take a buffer and create the missing groups localy void checkAndAddGroups(CommunicationBufferTemplated & buffer); /// register an elemental field to the given group name template inline dumper::Field * createElementalField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, const ElementTypeMap & nb_data_per_elem); /// register an elemental field to the given group name template inline dumper::Field * createElementalFilteredField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); #endif /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ElementGroups, element_groups, const ElementGroups &); const ElementGroup & getElementGroup(const std::string & name) const; const NodeGroup & getNodeGroup(const std::string & name) const; ElementGroup & getElementGroup(const std::string & name); NodeGroup & getNodeGroup(const std::string & name); UInt getNbElementGroups(UInt dimension = _all_dimensions) const; UInt getNbNodeGroups() { return node_groups.size(); }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id to create element and node groups ID id; /// memory_id to create element and node groups MemoryID memory_id; /// list of the node groups managed NodeGroups node_groups; /// list of the element groups managed ElementGroups element_groups; /// Mesh to which the element belongs const Mesh & mesh; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const GroupManager & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AKANTU_GROUP_MANAGER_HH__ */ diff --git a/src/model/dof_manager_default.cc b/src/model/dof_manager_default.cc index 40278f75e..617620eaa 100644 --- a/src/model/dof_manager_default.cc +++ b/src/model/dof_manager_default.cc @@ -1,920 +1,920 @@ /** * @file dof_manager_default.cc * * @author Nicolas Richart * * @date Tue Aug 11 16:21:01 2015 * * @brief Implementation of the default DOFManager * * @section LICENSE * * Copyright (©) 2010-2011 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 "dof_manager_default.hh" #include "communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include "node_synchronizer.hh" #include "non_linear_solver_default.hh" #include "sparse_matrix_aij.hh" #include "terms_to_assemble.hh" #include "time_step_solver_default.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = i; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { if (c_jcn >= c_irn) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addElementalMatrixToUnsymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(const ID & id, const MemoryID & memory_id) : DOFManager(id, memory_id), residual(0, 1, std::string(id + ":residual")), global_residual(nullptr), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), previous_global_blocked_dofs( 0, 1, std::string(id + ":previous_global_blocked_dofs")), dofs_type(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), global_equation_number(0, 1, "global_equation_number"), synchronizer(nullptr) {} /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(Mesh & mesh, const ID & id, const MemoryID & memory_id) : DOFManager(mesh, id, memory_id), residual(0, 1, std::string(id + ":residual")), global_residual(nullptr), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), previous_global_blocked_dofs( 0, 1, std::string(id + ":previous_global_blocked_dofs")), dofs_type(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), jacobian_release(0), global_equation_number(0, 1, "global_equation_number"), first_global_dof_id(0), synchronizer(nullptr) { if (this->mesh->isDistributed()) this->synchronizer = std::make_unique( *this, this->id + ":dof_synchronizer", this->memory_id); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::~DOFManagerDefault() = default; /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, Array & global_array, T scale_factor) { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = array_to_assemble.size() * array_to_assemble.getNbComponent(); AKANTU_DEBUG_ASSERT(equation_number.size() == nb_degree_of_freedoms, "The array to assemble does not have a correct size." << " (" << array_to_assemble.getID() << ")"); typename Array::const_scalar_iterator arr_it = array_to_assemble.begin_reinterpret(nb_degree_of_freedoms); Array::const_scalar_iterator equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++arr_it, ++equ_it) { global_array(*equ_it) += scale_factor * (*arr_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFDataDefault::DOFDataDefault(const ID & dof_id) : DOFData(dof_id), associated_nodes(0, 1, dof_id + "associated_nodes") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData & DOFManagerDefault::getNewDOFData(const ID & dof_id) { this->dofs[dof_id] = std::make_unique(dof_id); return *this->dofs[dof_id]; } /* -------------------------------------------------------------------------- */ class GlobalDOFInfoDataAccessor : public DataAccessor { public: using size_type = typename std::unordered_map>::size_type; GlobalDOFInfoDataAccessor() = default; void addDOFToNode(UInt node, UInt dof) { dofs_per_node[node].push_back(dof); } UInt getNthDOFForNode(UInt nth_dof, UInt node) const { return dofs_per_node.find(node)->second[nth_dof]; } UInt getNbData(const Array & nodes, const SynchronizationTag & tag) const override { if (tag == _gst_size) { return nodes.size() * sizeof(size_type); } if (tag == _gst_update) { UInt total_size = 0; for (auto node : nodes) { auto it = dofs_per_node.find(node); if (it != dofs_per_node.end()) total_size += CommunicationBuffer::sizeInBuffer(it->second); } return total_size; } return 0; } void packData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const override { for (auto node : nodes) { auto it = dofs_per_node.find(node); if (tag == _gst_size) { if (it != dofs_per_node.end()) { buffer << it->second.size(); } else { buffer << 0; } } else if (tag == _gst_update) { if (it != dofs_per_node.end()) buffer << it->second; } } } void unpackData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) override { for (auto node : nodes) { auto it = dofs_per_node.find(node); if (tag == _gst_size) { size_type size; buffer >> size; if (size != 0) dofs_per_node[node].resize(size); } else if (tag == _gst_update) { if (it != dofs_per_node.end()) buffer >> it->second; } } } protected: std::unordered_map> dofs_per_node; }; /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFsInternal(const ID & dof_id, UInt nb_dofs, UInt nb_pure_local_dofs) { // auto prank = this->communicator.whoAmI(); // auto psize = this->communicator.getNbProc(); // access the relevant data to update auto & dof_data = this->getDOFDataTyped(dof_id); const auto & support_type = dof_data.support_type; const auto & group = dof_data.group_support; if (support_type == _dst_nodal and group != "__mesh__") { auto & support_nodes = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); this->updateDOFsData( dof_data, nb_dofs, nb_pure_local_dofs, [&support_nodes](UInt node) -> UInt { return support_nodes[node]; }); } else { this->updateDOFsData(dof_data, nb_dofs, nb_pure_local_dofs, [](UInt node) -> UInt { return node; }); } // update the synchronizer if needed if (this->synchronizer) this->synchronizer->registerDOFs(dof_id); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { // stores the current numbers of dofs UInt nb_dofs_old = this->local_system_size; UInt nb_pure_local_dofs_old = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, support_type); UInt nb_dofs = this->local_system_size - nb_dofs_old; UInt nb_pure_local_dofs = this->pure_local_system_size - nb_pure_local_dofs_old; this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & group_support) { // stores the current numbers of dofs UInt nb_dofs_old = this->local_system_size; UInt nb_pure_local_dofs_old = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, group_support); UInt nb_dofs = this->local_system_size - nb_dofs_old; UInt nb_pure_local_dofs = this->pure_local_system_size - nb_pure_local_dofs_old; this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const MatrixType & matrix_type) { ID matrix_id = this->id + ":mtx:" + id; std::unique_ptr sm = std::make_unique(*this, matrix_type, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const ID & matrix_to_copy_id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatrixAIJ & sm_to_copy = this->getMatrix(matrix_to_copy_id); std::unique_ptr sm = std::make_unique(sm_to_copy, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } /* -------------------------------------------------------------------------- */ SparseMatrixAIJ & DOFManagerDefault::getMatrix(const ID & id) { SparseMatrix & matrix = DOFManager::getMatrix(id); return dynamic_cast(matrix); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManagerDefault::getNewNonLinearSolver(const ID & id, const NonLinearSolverType & type) { ID non_linear_solver_id = this->id + ":nls:" + id; std::unique_ptr nls; switch (type) { #if defined(AKANTU_IMPLICIT) case _nls_newton_raphson: case _nls_newton_raphson_modified: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } case _nls_linear: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } #endif case _nls_lumped: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } default: AKANTU_EXCEPTION("The asked type of non linear solver is not supported by " "this dof manager"); } return this->registerNonLinearSolver(non_linear_solver_id, nls); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerDefault::getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) { ID time_step_solver_id = this->id + ":tss:" + id; std::unique_ptr tss = std::make_unique( *this, type, non_linear_solver, time_step_solver_id, this->memory_id); return this->registerTimeStepSolver(time_step_solver_id, tss); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearResidual() { this->residual.resize(this->local_system_size); this->residual.clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearMatrix(const ID & mtx) { this->getMatrix(mtx).clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearLumpedMatrix(const ID & mtx) { this->getLumpedMatrix(mtx).clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateGlobalBlockedDofs() { auto it = this->dofs.begin(); auto end = this->dofs.end(); this->previous_global_blocked_dofs.copy(this->global_blocked_dofs); this->global_blocked_dofs.resize(this->local_system_size); this->global_blocked_dofs.clear(); for (; it != end; ++it) { if (!this->hasBlockedDOFs(it->first)) continue; DOFData & dof_data = *it->second; this->assembleToGlobalArray(it->first, *dof_data.blocked_dofs, this->global_blocked_dofs, true); } } /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::getArrayPerDOFs(const ID & dof_id, const Array & global_array, Array & local_array) const { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = equation_number.size(); local_array.resize(nb_degree_of_freedoms / local_array.getNbComponent()); auto loc_it = local_array.begin_reinterpret(nb_degree_of_freedoms); auto equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++loc_it, ++equ_it) { (*loc_it) = global_array(*equ_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void DOFManagerDefault::getEquationsNumbers(const ID & dof_id, // Array & equation_numbers) { // AKANTU_DEBUG_IN(); // this->getArrayPerDOFs(dof_id, this->global_equation_number, // equation_numbers); AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->global_solution, solution_array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getLumpedMatrix(lumped_mtx), lumped); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToResidual( const ID & dof_id, const Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); this->assembleToGlobalArray(dof_id, array_to_assemble, this->residual, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToLumpedMatrix( const ID & dof_id, const Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor) { AKANTU_DEBUG_IN(); Array & lumped = this->getLumpedMatrix(lumped_mtx); this->assembleToGlobalArray(dof_id, array_to_assemble, lumped, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { SparseMatrixAIJ & A = this->getMatrix(A_id); // Array data_cache(this->local_system_size, 1, 0.); this->data_cache.resize(this->local_system_size); this->data_cache.clear(); this->assembleToGlobalArray(dof_id, x, data_cache, 1.); A.matVecMul(data_cache, this->residual, scale_factor, 1.); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleLumpedMatMulVectToResidual( const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { const Array & A = this->getLumpedMatrix(A_id); // Array data_cache(this->local_system_size, 1, 0.); this->data_cache.resize(this->local_system_size); this->data_cache.clear(); this->assembleToGlobalArray(dof_id, x, data_cache, scale_factor); Array::const_scalar_iterator A_it = A.begin(); Array::const_scalar_iterator A_end = A.end(); Array::const_scalar_iterator x_it = data_cache.begin(); Array::scalar_iterator r_it = this->residual.begin(); for (; A_it != A_end; ++A_it, ++x_it, ++r_it) { *r_it += *A_it * *x_it; } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { AKANTU_DEBUG_IN(); DOFData & dof_data = this->getDOFData(dof_id); AKANTU_DEBUG_ASSERT(dof_data.support_type == _dst_nodal, "This function applies only on Nodal dofs"); this->addToProfile(matrix_id, dof_id, type, ghost_type); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt nb_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } else { if (dof_data.group_support != "__mesh__") { const Array & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); nb_element = group_elements.size(); filter_it = group_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } } AKANTU_DEBUG_ASSERT(elementary_mat.size() == nb_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good size."); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = dof_data.dof->getNbComponent(); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); auto conn_begin = connectivity.begin(nb_nodes_per_element); auto conn_it = conn_begin; UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; Vector element_eq_nb(nb_degree_of_freedom * nb_nodes_per_element); Array::const_matrix_iterator el_mat_it = elementary_mat.begin(size_mat, size_mat); for (UInt e = 0; e < nb_element; ++e, ++el_mat_it) { if (filter_it) conn_it = conn_begin + *filter_it; this->extractElementEquationNumber(equation_number, *conn_it, nb_degree_of_freedom, element_eq_nb); std::transform(element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (filter_it) ++filter_it; else ++conn_it; if (A.getMatrixType() == _symmetric) if (elemental_matrix_type == _symmetric) this->addSymmetricElementalMatrixToSymmetric(A, *el_mat_it, element_eq_nb, A.size()); else this->addUnsymmetricElementalMatrixToSymmetric(A, *el_mat_it, element_eq_nb, A.size()); else this->addElementalMatrixToUnsymmetric(A, *el_mat_it, element_eq_nb, A.size()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assemblePreassembledMatrix( const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) { const Array & equation_number_m = this->getLocalEquationNumbers(dof_id_m); const Array & equation_number_n = this->getLocalEquationNumbers(dof_id_n); SparseMatrixAIJ & A = this->getMatrix(matrix_id); for (const auto & term : terms) { UInt gi = this->localToGlobalEquationNumber(equation_number_m(term.i())); UInt gj = this->localToGlobalEquationNumber(equation_number_n(term.j())); A.add(gi, gj, term); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::addToProfile(const ID & matrix_id, const ID & dof_id, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const DOFData & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) return; auto mat_dof = std::make_pair(matrix_id, dof_id); auto type_pair = std::make_pair(type, ghost_type); auto prof_it = this->matrix_profiled_dofs.find(mat_dof); if (prof_it != this->matrix_profiled_dofs.end() && std::find(prof_it->second.begin(), prof_it->second.end(), type_pair) != prof_it->second.end()) return; UInt nb_degree_of_freedom_per_node = dof_data.dof->getNbComponent(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt size = A.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const auto & connectivity = this->mesh->getConnectivity(type, ghost_type); auto cbegin = connectivity.begin(nb_nodes_per_element); auto cit = cbegin; UInt nb_elements = connectivity.size(); UInt * ge_it = nullptr; if (dof_data.group_support != "__mesh__") { const Array & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); ge_it = group_elements.storage(); nb_elements = group_elements.size(); } UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom_per_node; Vector element_eq_nb(size_mat); for (UInt e = 0; e < nb_elements; ++e) { if (ge_it) cit = cbegin + *ge_it; this->extractElementEquationNumber( equation_number, *cit, nb_degree_of_freedom_per_node, element_eq_nb); std::transform(element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (ge_it) ++ge_it; else ++cit; for (UInt i = 0; i < size_mat; ++i) { UInt c_irn = element_eq_nb(i); if (c_irn < size) { for (UInt j = 0; j < size_mat; ++j) { UInt c_jcn = element_eq_nb(j); if (c_jcn < size) { A.add(c_irn, c_jcn); } } } } } this->matrix_profiled_dofs[mat_dof].push_back(type_pair); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::applyBoundary(const ID & matrix_id) { SparseMatrixAIJ & J = this->getMatrix(matrix_id); if (this->jacobian_release == J.getValueRelease()) { auto are_equal = std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), previous_global_blocked_dofs.begin()); - if (are_equal) + if (not are_equal) J.applyBoundary(); previous_global_blocked_dofs.copy(global_blocked_dofs); } else { J.applyBoundary(); } this->jacobian_release = J.getValueRelease(); } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getGlobalResidual() { if (this->synchronizer) { if (not this->global_residual) { this->global_residual = std::make_unique>(0, 1, "global_residual"); } if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->global_residual->resize(this->system_size); this->synchronizer->gather(this->residual, *this->global_residual); } else { this->synchronizer->gather(this->residual); } return *this->global_residual; } else { return this->residual; } } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getResidual() const { return this->residual; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::setGlobalSolution(const Array & solution) { if (this->synchronizer) { if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->synchronizer->scatter(this->global_solution, solution); } else { this->synchronizer->scatter(this->global_solution); } } else { AKANTU_DEBUG_ASSERT(solution.size() == this->global_solution.size(), "Sequential call to this function needs the solution " "to be the same size as the global_solution"); this->global_solution.copy(solution); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { DOFManager::onNodesAdded(nodes_list, event); if (this->synchronizer) this->synchronizer->onNodesAdded(nodes_list); } /* -------------------------------------------------------------------------- */ std::pair DOFManagerDefault::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { UInt nb_new_local_dofs, nb_new_pure_local; std::tie(nb_new_local_dofs, nb_new_pure_local) = DOFManager::updateNodalDOFs(dof_id, nodes_list); auto & dof_data = this->getDOFDataTyped(dof_id); updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local, [&nodes_list](UInt pos) -> UInt { return nodes_list[pos]; }); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateDOFsData( DOFDataDefault & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, const std::function & getNode) { auto prank = this->communicator.whoAmI(); auto psize = this->communicator.getNbProc(); // access the relevant data to update const auto & support_type = dof_data.support_type; GlobalDOFInfoDataAccessor data_accessor; // resize all relevant arrays this->residual.resize(this->local_system_size); this->dofs_type.resize(this->local_system_size); this->global_solution.resize(this->local_system_size); this->global_blocked_dofs.resize(this->local_system_size); this->previous_global_blocked_dofs.resize(this->local_system_size); this->global_equation_number.resize(this->local_system_size); for (auto & lumped_matrix : lumped_matrices) lumped_matrix.second->resize(this->local_system_size); dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_new_local_dofs); // determine the first local/global dof id to use Array nb_dofs_per_proc(psize); nb_dofs_per_proc(prank) = nb_new_pure_local; this->communicator.allGather(nb_dofs_per_proc); this->first_global_dof_id += std::accumulate( nb_dofs_per_proc.begin(), nb_dofs_per_proc.begin() + prank, 0); UInt first_dof_id = this->local_system_size - nb_new_local_dofs; if (support_type == _dst_nodal) { dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() + nb_new_local_dofs); } // update per dof info for (UInt d = 0; d < nb_new_local_dofs; ++d) { UInt local_eq_num = first_dof_id + d; dof_data.local_equation_number.push_back(local_eq_num); bool is_local_dof = true; // determine the dof type switch (support_type) { case _dst_nodal: { UInt node = getNode(d / dof_data.dof->getNbComponent()); this->dofs_type(local_eq_num) = this->mesh->getNodeType(node); dof_data.associated_nodes.push_back(node); is_local_dof = this->mesh->isLocalOrMasterNode(node); if (is_local_dof) { data_accessor.addDOFToNode(node, first_global_dof_id); } break; } case _dst_generic: { this->dofs_type(local_eq_num) = _nt_normal; break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } // update global id for local dofs if (is_local_dof) { this->global_equation_number(local_eq_num) = this->first_global_dof_id; this->global_to_local_mapping[this->first_global_dof_id] = local_eq_num; ++this->first_global_dof_id; } else { this->global_equation_number(local_eq_num) = 0; } } // synchronize the global numbering for slaves if (support_type == _dst_nodal && this->synchronizer) { auto nb_dofs_per_node = dof_data.dof->getNbComponent(); auto & node_synchronizer = this->mesh->getNodeSynchronizer(); node_synchronizer.synchronizeOnce(data_accessor, _gst_size); node_synchronizer.synchronizeOnce(data_accessor, _gst_update); std::vector counters(nb_new_local_dofs / nb_dofs_per_node); for (UInt d = 0; d < nb_new_local_dofs; ++d) { UInt local_eq_num = first_dof_id + d; if (not this->isSlaveDOF(local_eq_num)) continue; UInt node = d / nb_dofs_per_node; UInt dof_count = counters[node]++; node = getNode(node); UInt global_dof_id = data_accessor.getNthDOFForNode(dof_count, node); this->global_equation_number(local_eq_num) = global_dof_id; this->global_to_local_mapping[global_dof_id] = local_eq_num; } } } /* -------------------------------------------------------------------------- */ // register in factory static bool default_dof_manager_is_registered [[gnu::unused]] = DefaultDOFManagerFactory::getInstance().registerAllocator( "default", [](const ID & id, const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(id, mem_id); }); static bool dof_manager_is_registered [[gnu::unused]] = DOFManagerFactory::getInstance().registerAllocator( "default", [](Mesh & mesh, const ID & id, const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(mesh, id, mem_id); }); } // namespace akantu diff --git a/src/model/model.cc b/src/model/model.cc index 4b5cf6d36..7b0e963dd 100644 --- a/src/model/model.cc +++ b/src/model/model.cc @@ -1,367 +1,370 @@ /** * @file model.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Mon Oct 03 2011 * @date last modification: Thu Nov 19 2015 * * @brief implementation of model common parts * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "model.hh" #include "communicator.hh" #include "data_accessor.hh" #include "element_group.hh" #include "element_synchronizer.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Model::Model(Mesh & mesh, const ModelType & type, UInt dim, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), ModelSolver(mesh, type, id, memory_id), mesh(mesh), spatial_dimension(dim == _all_dimensions ? mesh.getSpatialDimension() : dim), is_pbc_slave_node(0, 1, "is_pbc_slave_node"), parser(getStaticParser()) { AKANTU_DEBUG_IN(); this->mesh.registerEventHandler(*this, _ehp_model); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Model::~Model() = default; /* -------------------------------------------------------------------------- */ // void Model::setParser(Parser & parser) { this->parser = &parser; } /* -------------------------------------------------------------------------- */ void Model::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); method = options.analysis_method; if (!this->hasDefaultSolver()) { this->initNewSolver(this->method); } initModel(); initFEEngineBoundary(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Model::initNewSolver(const AnalysisMethod & method) { ID solver_name; TimeStepSolverType tss_type; std::tie(solver_name, tss_type) = this->getDefaultSolverID(method); if (not this->hasSolver(solver_name)) { ModelSolverOptions options = this->getDefaultSolverOptions(tss_type); this->getNewSolver(solver_name, tss_type, options.non_linear_solver_type); for (auto && is_type : options.integration_scheme_type) { if (!this->hasIntegrationScheme(solver_name, is_type.first)) { this->setIntegrationScheme(solver_name, is_type.first, is_type.second, options.solution_type[is_type.first]); } } } this->method = method; this->setDefaultSolver(solver_name); } /* -------------------------------------------------------------------------- */ void Model::initPBC() { auto it = pbc_pair.begin(); auto end = pbc_pair.end(); is_pbc_slave_node.resize(mesh.getNbNodes()); #ifndef AKANTU_NDEBUG Array::const_vector_iterator coord_it = mesh.getNodes().begin(this->spatial_dimension); #endif while (it != end) { UInt i1 = (*it).first; is_pbc_slave_node(i1) = true; #ifndef AKANTU_NDEBUG UInt i2 = (*it).second; UInt slave = mesh.isDistributed() ? mesh.getGlobalNodesIds()(i1) : i1; UInt master = mesh.isDistributed() ? mesh.getGlobalNodesIds()(i2) : i2; AKANTU_DEBUG_INFO("pairing " << slave << " (" << Vector(coord_it[i1]) << ") with " << master << " (" << Vector(coord_it[i2]) << ")"); #endif ++it; } } /* -------------------------------------------------------------------------- */ void Model::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup(const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.dump(); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup(const std::string & group_name, const std::string & dumper_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup() { auto bit = mesh.element_group_begin(); auto bend = mesh.element_group_end(); for (; bit != bend; ++bit) { bit->second->dump(); } } /* -------------------------------------------------------------------------- */ void Model::setGroupDirectory(const std::string & directory) { auto bit = mesh.element_group_begin(); auto bend = mesh.element_group_end(); for (; bit != bend; ++bit) { bit->second->setDirectory(directory); } } /* -------------------------------------------------------------------------- */ void Model::setGroupDirectory(const std::string & directory, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.setDirectory(directory); } /* -------------------------------------------------------------------------- */ void Model::setGroupBaseName(const std::string & basename, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.setBaseName(basename); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Model::getGroupDumper(const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); return group.getDumper(); } /* -------------------------------------------------------------------------- */ // DUMPER stuff /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & field_id, dumper::Field * field, DumperIOHelper & dumper) { #ifdef AKANTU_USE_IOHELPER dumper.registerField(field_id, field); #endif } /* -------------------------------------------------------------------------- */ void Model::addDumpField(const std::string & field_id) { this->addDumpFieldToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldVector(const std::string & field_id) { this->addDumpFieldVectorToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldTensor(const std::string & field_id) { this->addDumpFieldTensorToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::setBaseName(const std::string & field_id) { mesh.setBaseName(field_id); } /* -------------------------------------------------------------------------- */ void Model::setBaseNameToDumper(const std::string & dumper_name, const std::string & basename) { mesh.setBaseNameToDumper(dumper_name, basename); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular, false); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id, group_name, _ek_regular, false); } /* -------------------------------------------------------------------------- */ void Model::removeDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void Model::removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.removeDumpFieldFromDumper(dumper_name, field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name, _ek_regular, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name, this->spatial_dimension, element_kind, padding_flag); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, UInt spatial_dimension, const ElementKind & element_kind, bool padding_flag) { #ifdef AKANTU_USE_IOHELPER dumper::Field * field = nullptr; if (!field) field = this->createNodalFieldReal(field_id, group_name, padding_flag); if (!field) field = this->createNodalFieldUInt(field_id, group_name, padding_flag); if (!field) field = this->createNodalFieldBool(field_id, group_name, padding_flag); if (!field) field = this->createElementalField(field_id, group_name, padding_flag, spatial_dimension, element_kind); if (!field) field = this->mesh.createFieldFromAttachedData(field_id, group_name, element_kind); if (!field) field = this->mesh.createFieldFromAttachedData(field_id, group_name, element_kind); - if (!field) +#ifndef AKANTU_NDEBUG + if (!field) { AKANTU_DEBUG_WARNING("No field could be found based on name: " << field_id); + } +#endif if (field) { DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name, group_name); this->addDumpGroupFieldToDumper(field_id, field, dumper); } #endif } /* -------------------------------------------------------------------------- */ void Model::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void Model::setDirectory(const std::string & directory) { mesh.setDirectory(directory); } /* -------------------------------------------------------------------------- */ void Model::setDirectoryToDumper(const std::string & dumper_name, const std::string & directory) { mesh.setDirectoryToDumper(dumper_name, directory); } /* -------------------------------------------------------------------------- */ void Model::setTextModeToDumper() { mesh.setTextModeToDumper(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/materials/internal_field.hh b/src/model/solid_mechanics/materials/internal_field.hh index 75d8254b2..d28b6a811 100644 --- a/src/model/solid_mechanics/materials/internal_field.hh +++ b/src/model/solid_mechanics/materials/internal_field.hh @@ -1,258 +1,258 @@ /** * @file internal_field.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Dec 08 2015 * * @brief Material internal properties * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_INTERNAL_FIELD_HH__ #define __AKANTU_INTERNAL_FIELD_HH__ namespace akantu { class Material; class FEEngine; /** * class for the internal fields of materials * to store values for each quadrature */ template class InternalField : public ElementTypeMapArray { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: InternalField(const ID & id, Material & material); ~InternalField() override; /// This constructor is only here to let cohesive elements compile InternalField(const ID & id, Material & material, FEEngine & fem, const ElementTypeMapArray & element_filter); /// More general constructor InternalField(const ID & id, Material & material, UInt dim, FEEngine & fem, const ElementTypeMapArray & element_filter); InternalField(const ID & id, const InternalField & other); private: InternalField operator=(__attribute__((unused)) const InternalField & other){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to reset the FEEngine for the internal field virtual void setFEEngine(FEEngine & fe_engine); /// function to reset the element kind for the internal virtual void setElementKind(ElementKind element_kind); /// initialize the field to a given number of component virtual void initialize(UInt nb_component); /// activate the history of this field virtual void initializeHistory(); /// resize the arrays and set the new element to 0 virtual void resize(); /// set the field to a given value v virtual void setDefaultValue(const T & v); /// reset all the fields to the default value virtual void reset(); /// save the current values in the history virtual void saveCurrentValues(); /// remove the quadrature points corresponding to suppressed elements virtual void removeIntegrationPoints(const ElementTypeMapArray & new_numbering); /// print the content - void printself(std::ostream & stream, int indent = 0) const override; + void printself(std::ostream & stream, int /*indent*/ = 0) const override; /// get the default value inline operator T() const; virtual FEEngine & getFEEngine() { return *fem; } virtual const FEEngine & getFEEngine() const { return *fem; } /// AKANTU_GET_MACRO(FEEngine, *fem, FEEngine &); protected: /// initialize the arrays in the ElementTypeMapArray void internalInitialize(UInt nb_component); /// set the values for new internals virtual void setArrayValues(T * begin, T * end); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: using type_iterator = typename ElementTypeMapArray::type_iterator; using filter_type_iterator = typename ElementTypeMapArray::type_iterator; /// get the type iterator on all types contained in the internal field type_iterator firstType(const GhostType & ghost_type = _not_ghost) const { return ElementTypeMapArray::firstType(this->spatial_dimension, ghost_type, this->element_kind); } /// get the type iterator on the last type contained in the internal field type_iterator lastType(const GhostType & ghost_type = _not_ghost) const { return ElementTypeMapArray::lastType(this->spatial_dimension, ghost_type, this->element_kind); } /// get the type iterator on all types contained in the internal field filter_type_iterator filterFirstType(const GhostType & ghost_type = _not_ghost) const { return this->element_filter.firstType(this->spatial_dimension, ghost_type, this->element_kind); } /// get the type iterator on the last type contained in the internal field filter_type_iterator filterLastType(const GhostType & ghost_type = _not_ghost) const { return this->element_filter.lastType(this->spatial_dimension, ghost_type, this->element_kind); } /// get the array for a given type of the element_filter const Array getFilter(const ElementType & type, const GhostType & ghost_type = _not_ghost) const { return this->element_filter(type, ghost_type); } /// get the Array corresponding to the type en ghost_type specified virtual Array & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost) { return ElementTypeMapArray::operator()(type, ghost_type); } virtual const Array & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost) const { return ElementTypeMapArray::operator()(type, ghost_type); } virtual Array & previous(const ElementType & type, const GhostType & ghost_type = _not_ghost) { AKANTU_DEBUG_ASSERT(previous_values != nullptr, "The history of the internal " << this->getID() << " has not been activated"); return this->previous_values->operator()(type, ghost_type); } virtual const Array & previous(const ElementType & type, const GhostType & ghost_type = _not_ghost) const { AKANTU_DEBUG_ASSERT(previous_values != nullptr, "The history of the internal " << this->getID() << " has not been activated"); return this->previous_values->operator()(type, ghost_type); } virtual InternalField & previous() { AKANTU_DEBUG_ASSERT(previous_values != nullptr, "The history of the internal " << this->getID() << " has not been activated"); return *(this->previous_values); } virtual const InternalField & previous() const { AKANTU_DEBUG_ASSERT(previous_values != nullptr, "The history of the internal " << this->getID() << " has not been activated"); return *(this->previous_values); } /// check if the history is used or not bool hasHistory() const { return (previous_values != nullptr); } /// get the kind treated by the internal const ElementKind & getElementKind() const { return element_kind; } /// return the number of components UInt getNbComponent() const { return nb_component; } /// return the spatial dimension corresponding to the internal element type /// loop filter UInt getSpatialDimension() const { return this->spatial_dimension; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the material for which this is an internal parameter Material & material; /// the fem containing the mesh and the element informations FEEngine * fem; /// Element filter if needed const ElementTypeMapArray & element_filter; /// default value T default_value; /// spatial dimension of the element to consider UInt spatial_dimension; /// ElementKind of the element to consider ElementKind element_kind; /// Number of component of the internal field UInt nb_component; /// Is the field initialized bool is_init; /// previous values InternalField * previous_values; }; /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const InternalField & _this) { _this.printself(stream); return stream; } } // akantu #endif /* __AKANTU_INTERNAL_FIELD_HH__ */ diff --git a/src/model/solid_mechanics/materials/internal_field_tmpl.hh b/src/model/solid_mechanics/materials/internal_field_tmpl.hh index 85147ea95..312fcaff2 100644 --- a/src/model/solid_mechanics/materials/internal_field_tmpl.hh +++ b/src/model/solid_mechanics/materials/internal_field_tmpl.hh @@ -1,320 +1,320 @@ /** * @file internal_field_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Tue Dec 08 2015 * * @brief Material internal properties * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_INTERNAL_FIELD_TMPL_HH__ #define __AKANTU_INTERNAL_FIELD_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ template InternalField::InternalField(const ID & id, Material & material) : ElementTypeMapArray(id, material.getID(), material.getMemoryID()), material(material), fem(&(material.getModel().getFEEngine())), element_filter(material.getElementFilter()), default_value(T()), spatial_dimension(material.getModel().getSpatialDimension()), element_kind(_ek_regular), nb_component(0), is_init(false), previous_values(nullptr) {} /* -------------------------------------------------------------------------- */ template InternalField::InternalField( const ID & id, Material & material, FEEngine & fem, const ElementTypeMapArray & element_filter) : ElementTypeMapArray(id, material.getID(), material.getMemoryID()), material(material), fem(&fem), element_filter(element_filter), default_value(T()), spatial_dimension(material.getSpatialDimension()), element_kind(_ek_regular), nb_component(0), is_init(false), previous_values(nullptr) {} /* -------------------------------------------------------------------------- */ template InternalField::InternalField( const ID & id, Material & material, UInt dim, FEEngine & fem, const ElementTypeMapArray & element_filter) : ElementTypeMapArray(id, material.getID(), material.getMemoryID()), material(material), fem(&fem), element_filter(element_filter), default_value(T()), spatial_dimension(dim), element_kind(_ek_regular), nb_component(0), is_init(false), previous_values(nullptr) {} /* -------------------------------------------------------------------------- */ template InternalField::InternalField(const ID & id, const InternalField & other) : ElementTypeMapArray(id, other.material.getID(), other.material.getMemoryID()), material(other.material), fem(other.fem), element_filter(other.element_filter), default_value(other.default_value), spatial_dimension(other.spatial_dimension), element_kind(other.element_kind), nb_component(other.nb_component), is_init(false), previous_values(nullptr) { AKANTU_DEBUG_ASSERT(other.is_init, "Cannot create a copy of a non initialized field"); this->internalInitialize(this->nb_component); } /* -------------------------------------------------------------------------- */ template InternalField::~InternalField() { if (this->is_init) { this->material.unregisterInternal(*this); } delete previous_values; } /* -------------------------------------------------------------------------- */ template void InternalField::setFEEngine(FEEngine & fe_engine) { this->fem = &fe_engine; } /* -------------------------------------------------------------------------- */ template void InternalField::setElementKind(ElementKind element_kind) { this->element_kind = element_kind; } /* -------------------------------------------------------------------------- */ template void InternalField::initialize(UInt nb_component) { internalInitialize(nb_component); } /* -------------------------------------------------------------------------- */ template void InternalField::initializeHistory() { if (!previous_values) previous_values = new InternalField("previous_" + this->getID(), *this); } /* -------------------------------------------------------------------------- */ template void InternalField::resize() { if (!this->is_init) return; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { filter_type_iterator it = this->filterFirstType(*gt); filter_type_iterator end = this->filterLastType(*gt); for (; it != end; ++it) { UInt nb_element = this->element_filter(*it, *gt).size(); UInt nb_quadrature_points = this->fem->getNbIntegrationPoints(*it, *gt); UInt new_size = nb_element * nb_quadrature_points; UInt old_size = 0; Array * vect = nullptr; if (this->exists(*it, *gt)) { vect = &(this->operator()(*it, *gt)); old_size = vect->size(); vect->resize(new_size); } else { vect = &(this->alloc(nb_element * nb_quadrature_points, nb_component, *it, *gt)); } this->setArrayValues(vect->storage() + old_size * vect->getNbComponent(), vect->storage() + new_size * vect->getNbComponent()); } } } /* -------------------------------------------------------------------------- */ template void InternalField::setDefaultValue(const T & value) { this->default_value = value; this->reset(); } /* -------------------------------------------------------------------------- */ template void InternalField::reset() { for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { type_iterator it = this->firstType(*gt); type_iterator end = this->lastType(*gt); for (; it != end; ++it) { Array & vect = this->operator()(*it, *gt); vect.clear(); this->setArrayValues(vect.storage(), vect.storage() + vect.size() * vect.getNbComponent()); } } } /* -------------------------------------------------------------------------- */ template void InternalField::internalInitialize(UInt nb_component) { if (!this->is_init) { this->nb_component = nb_component; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { filter_type_iterator it = this->filterFirstType(*gt); filter_type_iterator end = this->filterLastType(*gt); for (; it != end; ++it) { UInt nb_element = this->element_filter(*it, *gt).size(); UInt nb_quadrature_points = this->fem->getNbIntegrationPoints(*it, *gt); if (this->exists(*it, *gt)) this->operator()(*it, *gt).resize(nb_element * nb_quadrature_points); else this->alloc(nb_element * nb_quadrature_points, nb_component, *it, *gt); } } this->material.registerInternal(*this); this->is_init = true; } this->reset(); if (this->previous_values) this->previous_values->internalInitialize(nb_component); } /* -------------------------------------------------------------------------- */ template void InternalField::setArrayValues(T * begin, T * end) { for (; begin < end; ++begin) *begin = this->default_value; } /* -------------------------------------------------------------------------- */ template void InternalField::saveCurrentValues() { AKANTU_DEBUG_ASSERT(this->previous_values != nullptr, "The history of the internal " << this->getID() << " has not been activated"); if (!this->is_init) return; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { type_iterator it = this->firstType(*gt); type_iterator end = this->lastType(*gt); for (; it != end; ++it) { this->previous_values->operator()(*it, *gt) .copy(this->operator()(*it, *gt)); } } } /* -------------------------------------------------------------------------- */ template void InternalField::removeIntegrationPoints( const ElementTypeMapArray & new_numbering) { for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { ElementTypeMapArray::type_iterator it = new_numbering.firstType(_all_dimensions, *gt, _ek_not_defined); ElementTypeMapArray::type_iterator end = new_numbering.lastType(_all_dimensions, *gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if (this->exists(type, *gt)) { Array & vect = this->operator()(type, *gt); if (!vect.size()) continue; const Array & renumbering = new_numbering(type, *gt); UInt nb_quad_per_elem = fem->getNbIntegrationPoints(type, *gt); UInt nb_component = vect.getNbComponent(); Array tmp(renumbering.size() * nb_quad_per_elem, nb_component); AKANTU_DEBUG_ASSERT( tmp.size() == vect.size(), "Something strange append some mater was created from nowhere!!"); AKANTU_DEBUG_ASSERT( tmp.size() == vect.size(), "Something strange append some mater was created or disappeared in " << vect.getID() << "(" << vect.size() << "!=" << tmp.size() << ") " "!!"); UInt new_size = 0; for (UInt i = 0; i < renumbering.size(); ++i) { UInt new_i = renumbering(i); if (new_i != UInt(-1)) { memcpy(tmp.storage() + new_i * nb_component * nb_quad_per_elem, vect.storage() + i * nb_component * nb_quad_per_elem, nb_component * nb_quad_per_elem * sizeof(T)); ++new_size; } } tmp.resize(new_size * nb_quad_per_elem); vect.copy(tmp); } } } } /* -------------------------------------------------------------------------- */ template -void InternalField::printself(std::ostream & stream, int indent) const { +void InternalField::printself(std::ostream & stream, int indent [[gnu::unused]]) const { stream << "InternalField [ " << this->getID(); #if !defined(AKANTU_NDEBUG) if (AKANTU_DEBUG_TEST(dblDump)) { stream << std::endl; ElementTypeMapArray::printself(stream, indent + 3); } else { #endif stream << " {" << this->getData(_not_ghost).size() << " types - " << this->getData(_ghost).size() << " ghost types" << "}"; #if !defined(AKANTU_NDEBUG) } #endif stream << " ]"; } /* -------------------------------------------------------------------------- */ template <> inline void ParameterTyped >::setAuto( const ParserParameter & in_param) { Parameter::setAuto(in_param); Real r = in_param; param.setDefaultValue(r); } /* -------------------------------------------------------------------------- */ template inline InternalField::operator T() const { return default_value; } } // akantu #endif /* __AKANTU_INTERNAL_FIELD_TMPL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_python/material_python.cc b/src/model/solid_mechanics/materials/material_python/material_python.cc index 21367f690..24bd2dd46 100644 --- a/src/model/solid_mechanics/materials/material_python/material_python.cc +++ b/src/model/solid_mechanics/materials/material_python/material_python.cc @@ -1,164 +1,213 @@ /** * @file material_python.cc * * @author Guillaume Anciaux * * @date creation: Fri Nov 13 2015 * @date last modification: Fri Nov 13 2015 * * @brief Material python implementation * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_python.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MaterialPython::MaterialPython(SolidMechanicsModel & model, PyObject * obj, const ID & id) : Material(model, id), PythonFunctor(obj) { AKANTU_DEBUG_IN(); this->registerInternals(); std::vector param_names = - this->callFunctor >("registerParam"); + this->callFunctor>("registerParam"); for (UInt i = 0; i < param_names.size(); ++i) { std::stringstream sstr; sstr << "PythonParameter" << i; this->registerParam(param_names[i], local_params[param_names[i]], 0., - _pat_parsable, sstr.str()); + _pat_parsable | _pat_readable, sstr.str()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialPython::registerInternals() { std::vector internal_names = - this->callFunctor >("registerInternals"); + this->callFunctor>("registerInternals"); std::vector internal_sizes; try { internal_sizes = - this->callFunctor >("registerInternalSizes"); + this->callFunctor>("registerInternalSizes"); } catch (...) { internal_sizes.assign(internal_names.size(), 1); } for (UInt i = 0; i < internal_names.size(); ++i) { std::stringstream sstr; sstr << "PythonInternal" << i; this->internals[internal_names[i]] = - new InternalField(internal_names[i], *this); + std::make_unique>(internal_names[i], *this); AKANTU_DEBUG_INFO("alloc internal " << internal_names[i] << " " - << this->internals[internal_names[i]]); + << &this->internals[internal_names[i]]); this->internals[internal_names[i]]->initialize(internal_sizes[i]); } // making an internal with the quadrature points coordinates this->internals["quad_coordinates"] = - new InternalField("quad_coordinates", *this); + std::make_unique>("quad_coordinates", *this); auto && coords = *this->internals["quad_coordinates"]; coords.initialize(this->getSpatialDimension()); } /* -------------------------------------------------------------------------- */ + void MaterialPython::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); auto && coords = *this->internals["quad_coordinates"]; - this->model.getFEEngine().computeIntegrationPointsCoordinates(coords, - &this->element_filter); + this->model.getFEEngine().computeIntegrationPointsCoordinates( + coords, &this->element_filter); + + auto params = local_params; + params["rho"] = this->rho; + + try { + this->callFunctor("initMaterial", this->internals, params); + } catch (...) { + } + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ + void MaterialPython::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); + auto params = local_params; + params["rho"] = this->rho; + std::map *> internal_arrays; for (auto & i : this->internals) { auto & array = (*i.second)(el_type, ghost_type); auto & name = i.first; internal_arrays[name] = &array; } - auto params = local_params; - params["rho"] = this->rho; - this->callFunctor("computeStress", this->gradu(el_type, ghost_type), this->stress(el_type, ghost_type), internal_arrays, params); AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ -template -void MaterialPython::computeStress(Matrix & grad_u, Matrix & sigma, - std::vector & internal_iterators) { - std::vector inputs; - for (auto & i : internal_iterators) { - inputs.push_back(*i); - } - - for (UInt i = 0; i < inputs.size(); ++i) { - *internal_iterators[i] = inputs[i]; - } -} - /* -------------------------------------------------------------------------- */ void MaterialPython::computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { + auto params = local_params; + params["rho"] = this->rho; + std::map *> internal_arrays; for (auto & i : this->internals) { auto & array = (*i.second)(el_type, ghost_type); auto & name = i.first; internal_arrays[name] = &array; } - auto params = local_params; - params["rho"] = this->rho; - this->callFunctor("computeTangentModuli", this->gradu(el_type, ghost_type), tangent_matrix, internal_arrays, params); } /* -------------------------------------------------------------------------- */ Real MaterialPython::getPushWaveSpeed(const Element &) const { auto params = local_params; params["rho"] = this->rho; return this->callFunctor("getPushWaveSpeed", params); } +/* -------------------------------------------------------------------------- */ + +Real MaterialPython::getEnergyForType(const std::string & type, ElementType el_type) { + AKANTU_DEBUG_IN(); + + std::map *> internal_arrays; + for (auto & i : this->internals) { + auto & array = (*i.second)(el_type, _not_ghost); + auto & name = i.first; + internal_arrays[name] = &array; + } + + auto params = local_params; + params["rho"] = this->rho; + + auto & energy_density = *internal_arrays[type]; + + this->callFunctor("getEnergyDensity", type, energy_density, + this->gradu(el_type, _not_ghost), + this->stress(el_type, _not_ghost), internal_arrays, + params); + + Real energy = fem.integrate(energy_density, el_type, _not_ghost, + element_filter(el_type, _not_ghost)); + + AKANTU_DEBUG_OUT(); + return energy; +} + +/* -------------------------------------------------------------------------- */ + +Real MaterialPython::getEnergy(const std::string & type) { + AKANTU_DEBUG_IN(); + + if (this->internals.find(type) == this->internals.end()) { + AKANTU_EXCEPTION("unknown energy type: " + << type << " you must declare an internal named " << type); + } + + Real energy = 0.; + /// integrate the potential energy for each type of elements + Mesh::type_iterator it = element_filter.firstType(spatial_dimension); + Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension); + for (; it != last_type; ++it) { + energy += this->getEnergyForType(type, *it); + } + AKANTU_DEBUG_OUT(); + return energy; +} + +/* -------------------------------------------------------------------------- */ + } // akantu diff --git a/src/model/solid_mechanics/materials/material_python/material_python.hh b/src/model/solid_mechanics/materials/material_python/material_python.hh index 9237401dd..e675b338b 100644 --- a/src/model/solid_mechanics/materials/material_python/material_python.hh +++ b/src/model/solid_mechanics/materials/material_python/material_python.hh @@ -1,119 +1,97 @@ /** * @file material_python.hh * * @author Guillaume Anciaux * * @date creation: Fri Nov 13 2015 * @date last modification: Fri Nov 13 2015 * * @brief Python material * * @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 . * */ /** * @file material_python.hh * * @author Guillaume Anciaux * */ /* -------------------------------------------------------------------------- */ #include "python_functor.hh" /* -------------------------------------------------------------------------- */ #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_PYTHON_HH__ #define __AKANTU_MATERIAL_PYTHON_HH__ /* -------------------------------------------------------------------------- */ - - namespace akantu { class MaterialPython : public Material, PythonFunctor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialPython(SolidMechanicsModel & model, PyObject * obj, const ID & id = ""); ~MaterialPython() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void registerInternals(); void initMaterial() override; /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override; - /// constitutive law for a given quad point - template - void computeStress(Matrix & grad_u, Matrix & sigma, - std::vector & internal_iterators); - /// compute the tangent stiffness matrix for an element type void computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override; /// compute the push wave speed of the material Real getPushWaveSpeed(const Element & element) const override; -protected: - /// update the dissipated energy, must be called after the stress have been - /// computed - // virtual void updateEnergies(ElementType el_type, GhostType ghost_type){}; - - /// compute the tangent stiffness matrix for a given quadrature point - // inline void computeTangentModuliOnQuad(Matrix & tangent, Real & - // dam){}; + /// compute an energy of the material + Real getEnergy(const std::string & type) override; - /* ------------------------------------------------------------------------ */ - /* DataAccessor inherited members */ - /* ------------------------------------------------------------------------ */ -public: - /* ------------------------------------------------------------------------ */ - /* Accessors */ - /* ------------------------------------------------------------------------ */ -public: - // virtual Real getEnergy(std::string type){}; - // virtual Real getEnergy(std::string energy_id, ElementType type, UInt - // index){}; + /// compute an energy of the material + Real getEnergyForType(const std::string & type, ElementType el_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: std::map local_params; - std::map *> internals; + std::map>> internals; }; } // akantu #endif /* __AKANTU_MATERIAL_PYTHON_HH__ */ diff --git a/src/model/solid_mechanics/materials/random_internal_field_tmpl.hh b/src/model/solid_mechanics/materials/random_internal_field_tmpl.hh index a67c0a009..c5b3cb4bc 100644 --- a/src/model/solid_mechanics/materials/random_internal_field_tmpl.hh +++ b/src/model/solid_mechanics/materials/random_internal_field_tmpl.hh @@ -1,125 +1,125 @@ /** * @file random_internal_field_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Tue Dec 08 2015 * * @brief Random internal material parameter implementation * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_random_generator.hh" #include "internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_RANDOM_INTERNAL_FIELD_TMPL_HH__ #define __AKANTU_RANDOM_INTERNAL_FIELD_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ template class BaseField, template class Generator> RandomInternalField::RandomInternalField( const ID & id, Material & material) : BaseField(id, material), random_parameter(T()) {} /* -------------------------------------------------------------------------- */ template class BaseField, template class Generator> RandomInternalField::~RandomInternalField() = default; /* -------------------------------------------------------------------------- */ template class BaseField, template class Generator> void RandomInternalField::initialize( UInt nb_component) { this->internalInitialize(nb_component); } /* ------------------------------------------------------------------------ */ template class BaseField, template class Generator> void RandomInternalField::setDefaultValue( const T & value) { random_parameter.setBaseValue(value); this->reset(); } /* ------------------------------------------------------------------------ */ template class BaseField, template class Generator> void RandomInternalField::setRandomDistribution( const RandomParameter & param) { random_parameter = param; this->reset(); } /* ------------------------------------------------------------------------ */ template class BaseField, template class Generator> void RandomInternalField::printself( - std::ostream & stream, int indent) const { + std::ostream & stream, int indent [[gnu::unused]]) const { stream << "RandomInternalField [ "; random_parameter.printself(stream); stream << " ]"; #if !defined(AKANTU_NDEBUG) if (AKANTU_DEBUG_TEST(dblDump)) { stream << std::endl; InternalField::printself(stream, indent); } #endif } /* -------------------------------------------------------------------------- */ template class BaseField, template class Generator> void RandomInternalField::setArrayValues(T * begin, T * end) { random_parameter.template setValues(begin, end); } /* -------------------------------------------------------------------------- */ template class BaseField, template class Generator> inline RandomInternalField::operator Real() const { return random_parameter.getBaseValue(); } /* -------------------------------------------------------------------------- */ template <> inline void ParameterTyped>::setAuto( const ParserParameter & in_param) { Parameter::setAuto(in_param); RandomParameter r = in_param; param.setRandomDistribution(r); } /* -------------------------------------------------------------------------- */ } // akantu #endif /* __AKANTU_RANDOM_INTERNAL_FIELD_TMPL_HH__ */ diff --git a/src/python/python_functor.cc b/src/python/python_functor.cc index bdfd9d132..15986dcf1 100644 --- a/src/python/python_functor.cc +++ b/src/python/python_functor.cc @@ -1,79 +1,80 @@ /** * @file python_functor.cc * * @author Guillaume Anciaux * * @date creation: Thu Feb 21 2013 * @date last modification: Fri Nov 13 2015 * * @brief Python functor interface * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "python_functor.hh" #include "aka_common.hh" +#include "internal_field.hh" #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ PythonFunctor::PythonFunctor(PyObject * obj) : python_obj(obj) {} /* -------------------------------------------------------------------------- */ PyObject * PythonFunctor::callFunctor(PyObject * functor, PyObject * args, PyObject * kwargs) const { if (!PyCallable_Check(functor)) AKANTU_EXCEPTION("Provided functor is not a function"); PyObject * pValue = PyObject_Call(functor, args, kwargs); PyObject * exception_type = PyErr_Occurred(); if (exception_type) { PyObject * exception; PyObject * traceback; PyErr_Fetch(&exception_type, &exception, &traceback); PyObject_Print(exception_type, stdout, Py_PRINT_RAW); PyObject_Print(exception, stdout, Py_PRINT_RAW); std::stringstream sstr; sstr << "Exception occured while calling the functor: "; PyObject * exception_mesg = PyObject_GetAttrString(exception, "message"); #if PY_MAJOR_VERSION >= 3 if (exception_mesg && PyUnicode_Check(exception_mesg)) #else if (exception_mesg && PyString_Check(exception_mesg)) #endif sstr << this->convertToAkantu(exception_mesg); else sstr << this->convertToAkantu(exception); AKANTU_EXCEPTION(sstr.str()); } return pValue; } /* -------------------------------------------------------------------------- */ } // akantu diff --git a/src/python/python_functor.hh b/src/python/python_functor.hh index 7516baeea..22bc4a2c5 100644 --- a/src/python/python_functor.hh +++ b/src/python/python_functor.hh @@ -1,129 +1,143 @@ /** * @file python_functor.hh * * @author Guillaume Anciaux * * @date creation: Fri Jun 18 2010 * @date last modification: Sun Nov 15 2015 * * @brief Python Functor interface * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" +#include "internal_field.hh" +#include "element_group.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PYTHON_FUNCTOR_HH__ #define __AKANTU_PYTHON_FUNCTOR_HH__ namespace akantu { class PythonFunctor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: PythonFunctor(PyObject * obj); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ -protected: + /// call the python functor PyObject * callFunctor(PyObject * functor, PyObject * args, PyObject * kwargs) const; /// call the python functor from variadic types template return_type callFunctor(const std::string & functor_name, Params &... parameters) const; /// empty function to cose the recursive template loop - inline void packArguments(std::vector & pArgs) const; + static inline void packArguments(std::vector & pArgs); /// get the python function object inline PyObject * getPythonFunction(const std::string & functor_name) const; /// variadic template for unknown number of arguments to unpack template - inline void packArguments(std::vector & pArgs, T & p, - Args &... params) const; + static inline void packArguments(std::vector & pArgs, T & p, + Args &... params); /// convert an akantu object to python template - inline PyObject * convertToPython(const T & akantu_obj) const; + static inline PyObject * convertToPython(const T & akantu_obj); /// convert a stl vector to python template - inline PyObject * convertToPython(const std::vector & akantu_obj) const; + static inline PyObject * convertToPython(const std::vector & akantu_obj); /// convert a stl vector to python template - inline PyObject * convertToPython(const std::vector *> & akantu_obj) const; + static inline PyObject * convertToPython(const std::vector *> & akantu_obj); + /// convert a stl vector to python + template + static inline PyObject * convertToPython(const std::unique_ptr & akantu_obj); + /// convert a stl vector to python template - inline PyObject * convertToPython(const std::map & akantu_obj) const; + static inline PyObject * convertToPython(const std::map & akantu_obj); /// convert an akantu vector to python template - inline PyObject * convertToPython(const Vector & akantu_obj) const; + static inline PyObject * convertToPython(const Vector & akantu_obj); + /// convert an akantu internal to python + template + static inline PyObject * convertToPython(const InternalField & akantu_obj); + + /// convert an akantu vector to python + template + static inline PyObject * convertToPython(const ElementTypeMapArray & akantu_obj); + /// convert an akantu vector to python template - inline PyObject * convertToPython(const Array & akantu_obj) const; + static inline PyObject * convertToPython(const Array & akantu_obj); /// convert an akantu vector to python template - inline PyObject * convertToPython(Array * akantu_obj) const; + static inline PyObject * convertToPython(Array * akantu_obj); /// convert a akantu matrix to python template - inline PyObject * convertToPython(const Matrix & akantu_obj) const; + static inline PyObject * convertToPython(const Matrix & akantu_obj); /// convert a python object to an akantu object template - inline return_type convertToAkantu(PyObject * python_obj) const; + static inline return_type convertToAkantu(PyObject * python_obj); /// convert a python object to an akantu object template - inline std::vector convertListToAkantu(PyObject * python_obj) const; + static inline std::vector convertListToAkantu(PyObject * python_obj); /// returns the numpy data type code - template inline int getPythonDataTypeCode() const; + template static inline int getPythonDataTypeCode(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: PyObject * python_obj; }; } // akantu #endif /* __AKANTU_PYTHON_FUNCTOR_HH__ */ #include "python_functor_inline_impl.cc" diff --git a/src/python/python_functor_inline_impl.cc b/src/python/python_functor_inline_impl.cc index 88ed09184..38ceb9b6c 100644 --- a/src/python/python_functor_inline_impl.cc +++ b/src/python/python_functor_inline_impl.cc @@ -1,351 +1,463 @@ /** * @file python_functor_inline_impl.cc * * @author Guillaume Anciaux * * @date creation: Fri Nov 13 2015 * @date last modification: Wed Nov 18 2015 * * @brief Python functor interface * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "integration_point.hh" +#include "internal_field.hh" /* -------------------------------------------------------------------------- */ #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include #include #include #if PY_MAJOR_VERSION >= 3 #include #endif /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PYTHON_FUNCTOR_INLINE_IMPL_CC__ #define __AKANTU_PYTHON_FUNCTOR_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ -template inline int PythonFunctor::getPythonDataTypeCode() const { +template inline int PythonFunctor::getPythonDataTypeCode() { AKANTU_EXCEPTION("undefined type: " << debug::demangle(typeid(T).name())); } /* -------------------------------------------------------------------------- */ -template <> inline int PythonFunctor::getPythonDataTypeCode() const { +template <> inline int PythonFunctor::getPythonDataTypeCode() { int data_typecode = NPY_NOTYPE; size_t s = sizeof(bool); switch (s) { case 1: data_typecode = NPY_BOOL; break; case 2: data_typecode = NPY_UINT16; break; case 4: data_typecode = NPY_UINT32; break; case 8: data_typecode = NPY_UINT64; break; } return data_typecode; } /* -------------------------------------------------------------------------- */ -template <> inline int PythonFunctor::getPythonDataTypeCode() const { +template <> inline int PythonFunctor::getPythonDataTypeCode() { return NPY_DOUBLE; } /* -------------------------------------------------------------------------- */ -template -PyObject * PythonFunctor::convertToPython(const T &) const { - AKANTU_DEBUG_ERROR(__func__ << " : not implemented yet !" << std::endl - << debug::demangle(typeid(T).name())); +template <> inline int PythonFunctor::getPythonDataTypeCode() { + return NPY_UINT; } + /* -------------------------------------------------------------------------- */ template <> inline PyObject * -PythonFunctor::convertToPython(const double & akantu_object) const { +PythonFunctor::convertToPython(const double & akantu_object) { return PyFloat_FromDouble(akantu_object); } /* -------------------------------------------------------------------------- */ template <> inline PyObject * -PythonFunctor::convertToPython(const UInt & akantu_object) const { +PythonFunctor::convertToPython(const UInt & akantu_object) { #if PY_MAJOR_VERSION >= 3 return PyLong_FromLong(akantu_object); #else return PyInt_FromLong(akantu_object); #endif } /* -------------------------------------------------------------------------- */ + template <> inline PyObject * -PythonFunctor::convertToPython(const bool & akantu_object) const { - return PyBool_FromLong(long(akantu_object)); +PythonFunctor::convertToPython(const int & akantu_object) { +#if PY_MAJOR_VERSION >= 3 + return PyLong_FromLong(akantu_object); +#else + return PyInt_FromLong(akantu_object); +#endif } /* -------------------------------------------------------------------------- */ + +template <> +inline PyObject * +PythonFunctor::convertToPython(const bool & akantu_object) { + return PyBool_FromLong(long(akantu_object)); +} +/* -------------------------------------------------------------------------- */ + +template <> +inline PyObject * +PythonFunctor::convertToPython(const std::string & str) { +#if PY_MAJOR_VERSION >= 3 + return PyUnicode_FromString(str.c_str()); +#else + return PyString_FromString(str.c_str()); +#endif +} + +/* -------------------------------------------------------------------------- + */ + +template <> +inline PyObject * +PythonFunctor::convertToPython(const NodeGroup & group) { + return PythonFunctor::convertToPython(group.getNodes()); +} + +/* -------------------------------------------------------------------------- + */ + +template <> +inline PyObject * +PythonFunctor::convertToPython(const ElementGroup & group) { + + PyObject * res = PyDict_New(); + + PyDict_SetItem(res, + PythonFunctor::convertToPython(std::string("element_group")), + PythonFunctor::convertToPython(group.getElements())); + + PyDict_SetItem(res, PythonFunctor::convertToPython(std::string("node_group")), + PythonFunctor::convertToPython(group.getNodeGroup())); + + return res; +} + +/* -------------------------------------------------------------------------- + */ + +template <> +inline PyObject * +PythonFunctor::convertToPython(ElementGroup * const & group) { + return PythonFunctor::convertToPython(*group); +} + +/* -------------------------------------------------------------------------- + */ + +template <> +inline PyObject * +PythonFunctor::convertToPython(const ElementType & type) { + // std::stringstream sstr; + // sstr << type; + // return PythonFunctor::convertToPython(sstr.str()); + return PythonFunctor::convertToPython(int(type)); +} + +/* -------------------------------------------------------------------------- + */ + template inline PyObject * -PythonFunctor::convertToPython(const std::vector & array) const { +PythonFunctor::convertToPython(const ElementTypeMapArray & map) { + + std::map *> res; + for (const auto & type : map.elementTypes()) { + std::stringstream sstr; + sstr << type; + res[sstr.str()] = const_cast *>(&map(type)); + } + return PythonFunctor::convertToPython(res); +} + +/* -------------------------------------------------------------------------- + */ + +template PyObject * PythonFunctor::convertToPython(const T &) { + AKANTU_EXCEPTION(__PRETTY_FUNCTION__ << " : not implemented yet !" + << std::endl + << debug::demangle(typeid(T).name())); +} + +/* -------------------------------------------------------------------------- */ +template +inline PyObject * PythonFunctor::convertToPython(const std::vector & array) { int data_typecode = getPythonDataTypeCode(); npy_intp dims[1] = {int(array.size())}; PyObject * obj = PyArray_SimpleNewFromData(1, dims, data_typecode, const_cast(&array[0])); auto * res = (PyArrayObject *)obj; return (PyObject *)res; } /* -------------------------------------------------------------------------- */ template inline PyObject * -PythonFunctor::convertToPython(const std::vector *> & array) const { +PythonFunctor::convertToPython(const std::vector *> & array) { PyObject * res = PyDict_New(); for (auto a : array) { - PyObject * obj = this->convertToPython(*a); - PyObject * name = this->convertToPython(a->getID()); + PyObject * obj = PythonFunctor::convertToPython(*a); + PyObject * name = PythonFunctor::convertToPython(a->getID()); PyDict_SetItem(res, name, obj); } return (PyObject *)res; } /* -------------------------------------------------------------------------- */ template -inline PyObject * -PythonFunctor::convertToPython(const std::map & map) const { +inline PyObject * PythonFunctor::convertToPython(const std::map & map) { PyObject * res = PyDict_New(); - for (auto a : map) { - PyObject * key = this->convertToPython(a.first); - PyObject * value = this->convertToPython(a.second); + for (auto && a : map) { + PyObject * key = PythonFunctor::convertToPython(a.first); + PyObject * value = PythonFunctor::convertToPython(a.second); PyDict_SetItem(res, key, value); } return (PyObject *)res; } /* -------------------------------------------------------------------------- */ template -PyObject * PythonFunctor::convertToPython(const Vector & array) const { +PyObject * PythonFunctor::convertToPython(const Vector & array) { int data_typecode = getPythonDataTypeCode(); npy_intp dims[1] = {array.size()}; PyObject * obj = PyArray_SimpleNewFromData(1, dims, data_typecode, array.storage()); auto * res = (PyArrayObject *)obj; return (PyObject *)res; } /* -------------------------------------------------------------------------- */ template -PyObject * PythonFunctor::convertToPython(const Array & array) const { +inline PyObject * +PythonFunctor::convertToPython(const InternalField & internals) { + return convertToPython( + static_cast &>(internals)); +} + +/* --------------------------------------------------------------------- */ + +template +inline PyObject * PythonFunctor::convertToPython(const std::unique_ptr & u) { + return convertToPython(*u); +} + +/* --------------------------------------------------------------------- */ + +template +PyObject * PythonFunctor::convertToPython(const Array & array) { int data_typecode = getPythonDataTypeCode(); npy_intp dims[2] = {array.size(), array.getNbComponent()}; PyObject * obj = PyArray_SimpleNewFromData(2, dims, data_typecode, array.storage()); auto * res = (PyArrayObject *)obj; return (PyObject *)res; } -/* -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ template -PyObject * PythonFunctor::convertToPython(Array * array) const { - return this->convertToPython(*array); +PyObject * PythonFunctor::convertToPython(Array * array) { + return PythonFunctor::convertToPython(*array); } -/* -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ template -PyObject * PythonFunctor::convertToPython(const Matrix & mat) const { +PyObject * PythonFunctor::convertToPython(const Matrix & mat) { int data_typecode = getPythonDataTypeCode(); npy_intp dims[2] = {mat.size(0), mat.size(1)}; PyObject * obj = PyArray_SimpleNewFromData(2, dims, data_typecode, mat.storage()); auto * res = (PyArrayObject *)obj; return (PyObject *)res; } -/* -------------------------------------------------------------------------- */ -template <> -inline PyObject * -PythonFunctor::convertToPython(const std::string & str) const { -#if PY_MAJOR_VERSION >= 3 - return PyUnicode_FromString(str.c_str()); -#else - return PyString_FromString(str.c_str()); -#endif -} +/* ---------------------------------------------------------------------- */ -/* -------------------------------------------------------------------------- */ template <> -inline PyObject * PythonFunctor::convertToPython( - const IntegrationPoint & qp) const { +inline PyObject * +PythonFunctor::convertToPython(const IntegrationPoint & qp) { PyObject * input = PyDict_New(); - PyObject * num_point = this->convertToPython(qp.num_point); - PyObject * global_num = this->convertToPython(qp.global_num); - PyObject * material_id = this->convertToPython(qp.material_id); - PyObject * position = this->convertToPython(qp.getPosition()); + PyObject * num_point = PythonFunctor::convertToPython(qp.num_point); + PyObject * global_num = PythonFunctor::convertToPython(qp.global_num); + PyObject * material_id = PythonFunctor::convertToPython(qp.material_id); + PyObject * position = PythonFunctor::convertToPython(qp.getPosition()); PyDict_SetItemString(input, "num_point", num_point); PyDict_SetItemString(input, "global_num", global_num); PyDict_SetItemString(input, "material_id", material_id); PyDict_SetItemString(input, "position", position); return input; } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ inline PyObject * PythonFunctor::getPythonFunction(const std::string & functor_name) const { #if PY_MAJOR_VERSION < 3 if (!PyInstance_Check(this->python_obj)) AKANTU_EXCEPTION("Python object is not an instance"); #else // does not make sense to check everything is an instance of object in python 3 #endif if (not PyObject_HasAttrString(this->python_obj, functor_name.c_str())) AKANTU_EXCEPTION("Python dictionary has no " << functor_name << " entry"); PyObject * pFunctor = PyObject_GetAttrString(this->python_obj, functor_name.c_str()); return pFunctor; } -/* -------------------------------------------------------------------------- */ -inline void -PythonFunctor::packArguments(__attribute__((unused)) - std::vector & p_args) const {} +/* -------------------------------------------------------------------------- + */ +inline void PythonFunctor::packArguments(__attribute__((unused)) + std::vector & p_args) {} -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ template inline void PythonFunctor::packArguments(std::vector & p_args, - T & p, Args &... params) const { - p_args.push_back(this->convertToPython(p)); + T & p, Args &... params) { + p_args.push_back(PythonFunctor::convertToPython(p)); if (sizeof...(params) != 0) - this->packArguments(p_args, params...); + PythonFunctor::packArguments(p_args, params...); } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ template return_type PythonFunctor::callFunctor(const std::string & functor_name, Params &... parameters) const { _import_array(); std::vector arg_vector; this->packArguments(arg_vector, parameters...); PyObject * pArgs = PyTuple_New(arg_vector.size()); for (UInt i = 0; i < arg_vector.size(); ++i) { PyTuple_SetItem(pArgs, i, arg_vector[i]); } PyObject * kwargs = PyDict_New(); - PyObject * pFunctor = getPythonFunction(functor_name); + PyObject * pFunctor = this->getPythonFunction(functor_name); PyObject * res = this->callFunctor(pFunctor, pArgs, kwargs); Py_XDECREF(pArgs); Py_XDECREF(kwargs); return this->convertToAkantu(res); } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ template -inline return_type PythonFunctor::convertToAkantu(PyObject * python_obj) const { +inline return_type PythonFunctor::convertToAkantu(PyObject * python_obj) { if (PyList_Check(python_obj)) { - return this->convertListToAkantu( + return PythonFunctor::convertListToAkantu( python_obj); } AKANTU_DEBUG_TO_IMPLEMENT(); } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ template <> -inline void PythonFunctor::convertToAkantu(PyObject * python_obj) const { +inline void PythonFunctor::convertToAkantu(PyObject * python_obj) { if (python_obj != Py_None) AKANTU_DEBUG_WARNING( "functor return a value while none was expected: ignored"); } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ template <> inline std::string -PythonFunctor::convertToAkantu(PyObject * python_obj) const { +PythonFunctor::convertToAkantu(PyObject * python_obj) { #if PY_MAJOR_VERSION >= 3 if (!PyUnicode_Check(python_obj)) AKANTU_EXCEPTION("cannot convert object to string"); std::wstring unicode_str(PyUnicode_AsWideCharString(python_obj, NULL)); std::wstring_convert, wchar_t> converter; return converter.to_bytes(unicode_str); #else if (!PyString_Check(python_obj)) AKANTU_EXCEPTION("cannot convert object to string"); return PyString_AsString(python_obj); #endif } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ template <> -inline Real PythonFunctor::convertToAkantu(PyObject * python_obj) const { +inline Real PythonFunctor::convertToAkantu(PyObject * python_obj) { if (!PyFloat_Check(python_obj)) AKANTU_EXCEPTION("cannot convert object to float"); return PyFloat_AsDouble(python_obj); } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ template <> -inline UInt PythonFunctor::convertToAkantu(PyObject * python_obj) const { +inline UInt PythonFunctor::convertToAkantu(PyObject * python_obj) { #if PY_MAJOR_VERSION >= 3 if (!PyLong_Check(python_obj)) AKANTU_EXCEPTION("cannot convert object to integer"); return PyLong_AsLong(python_obj); #else if (!PyInt_Check(python_obj)) AKANTU_EXCEPTION("cannot convert object to integer"); return PyInt_AsLong(python_obj); #endif } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ template inline std::vector -PythonFunctor::convertListToAkantu(PyObject * python_obj) const { +PythonFunctor::convertListToAkantu(PyObject * python_obj) { std::vector res; UInt size = PyList_Size(python_obj); for (UInt i = 0; i < size; ++i) { PyObject * item = PyList_GET_ITEM(python_obj, i); - res.push_back(this->convertToAkantu(item)); + res.push_back(PythonFunctor::convertToAkantu(item)); } return res; } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ } // akantu #endif /* __AKANTU_PYTHON_FUNCTOR_INLINE_IMPL_CC__ */ diff --git a/src/synchronizer/element_synchronizer.cc b/src/synchronizer/element_synchronizer.cc index 25a7f77d1..9ffc5f5f8 100644 --- a/src/synchronizer/element_synchronizer.cc +++ b/src/synchronizer/element_synchronizer.cc @@ -1,349 +1,353 @@ /** * @file element_synchronizer.cc * * @author Guillaume Anciaux * @author Dana Christen * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Sep 01 2010 * @date last modification: Fri Jan 22 2016 * * @brief implementation of a communicator using a static_communicator for * real * send/receive * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "element_synchronizer.hh" #include "aka_common.hh" #include "mesh.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ElementSynchronizer::ElementSynchronizer(Mesh & mesh, const ID & id, MemoryID memory_id, bool register_to_event_manager, EventHandlerPriority event_priority) : SynchronizerImpl(mesh.getCommunicator(), id, memory_id), mesh(mesh), element_to_prank("element_to_prank", id, memory_id) { AKANTU_DEBUG_IN(); if (register_to_event_manager) this->mesh.registerEventHandler(*this, event_priority); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementSynchronizer::~ElementSynchronizer() = default; /* -------------------------------------------------------------------------- */ void ElementSynchronizer::substituteElements( const std::map & old_to_new_elements) { auto found_element_end = old_to_new_elements.end(); // substitute old elements with new ones for (auto && sr : iterate_send_recv) { for (auto && scheme_pair : communications.iterateSchemes(sr)) { auto & list = scheme_pair.second; for (auto & el : list) { auto found_element_it = old_to_new_elements.find(el); if (found_element_it != found_element_end) el = found_element_it->second; } } } } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::onElementsChanged( const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray &, const ChangedElementsEvent &) { // create a map to link old elements to new ones std::map old_to_new_elements; for (UInt el = 0; el < old_elements_list.size(); ++el) { AKANTU_DEBUG_ASSERT(old_to_new_elements.find(old_elements_list(el)) == old_to_new_elements.end(), "The same element cannot appear twice in the list"); old_to_new_elements[old_elements_list(el)] = new_elements_list(el); } substituteElements(old_to_new_elements); communications.invalidateSizes(); } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::onElementsRemoved( const Array & element_to_remove, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent &) { AKANTU_DEBUG_IN(); this->removeElements(element_to_remove); this->renumberElements(new_numbering); communications.invalidateSizes(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::buildElementToPrank() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); element_to_prank.initialize(mesh, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined, _with_nb_element = true, _default_value = rank); /// assign prank to all ghost elements for (auto && scheme : communications.iterateSchemes(_recv)) { auto & recv = scheme.second; auto proc = scheme.first; for (auto & element : recv) { element_to_prank(element) = proc; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Int ElementSynchronizer::getRank(const Element & element) const { if (not element_to_prank.exists(element.type, element.ghost_type)) { // Nicolas: Ok This is nasty I know.... const_cast(this)->buildElementToPrank(); } return element_to_prank(element); } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::reset() { AKANTU_DEBUG_IN(); communications.resetSchemes(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::removeElements( const Array & element_to_remove) { AKANTU_DEBUG_IN(); std::vector send_requests; std::map> list_of_elements_per_proc; auto filter_list = [](const Array & keep, Array & list) { Array new_list; for (UInt e = 0; e < keep.size() - 1; ++e) { Element & el = list(keep(e)); new_list.push_back(el); } list.copy(new_list); }; // Handling ghost elements for (auto && scheme_pair : communications.iterateSchemes(_recv)) { auto proc = scheme_pair.first; auto & recv = scheme_pair.second; Array & keep_list = list_of_elements_per_proc[proc]; auto rem_it = element_to_remove.begin(); auto rem_end = element_to_remove.end(); for (auto && pair : enumerate(recv)) { const auto & el = std::get<1>(pair); auto pos = std::find(rem_it, rem_end, el); if (pos == rem_end) { keep_list.push_back(std::get<0>(pair)); } } keep_list.push_back(UInt(-1)); // To no send empty arrays auto && tag = Tag::genTag(proc, 0, Tag::_MODIFY_SCHEME, this->hash_id); AKANTU_DEBUG_INFO("Sending a message of size " << keep_list.size() << " to proc " << proc << " TAG(" << tag << ")"); send_requests.push_back(this->communicator.asyncSend(keep_list, proc, tag)); +#ifndef AKANTU_NDEBUG auto old_size = recv.size(); +#endif filter_list(keep_list, recv); AKANTU_DEBUG_INFO("I had " << old_size << " elements to recv from proc " << proc << " and " << keep_list.size() << " elements to keep. I have " << recv.size() << " elements left."); } for (auto && scheme_pair : communications.iterateSchemes(_send)) { auto proc = scheme_pair.first; auto & send = scheme_pair.second; CommunicationStatus status; auto && tag = Tag::genTag(rank, 0, Tag::_MODIFY_SCHEME, this->hash_id); AKANTU_DEBUG_INFO("Getting number of elements of proc " << proc << " to keep - TAG(" << tag << ")"); this->communicator.probe(proc, tag, status); Array keep_list(status.size()); AKANTU_DEBUG_INFO("Receiving list of elements (" << keep_list.size() << " elements) to keep for proc " << proc << " TAG(" << tag << ")"); this->communicator.receive(keep_list, proc, tag); +#ifndef AKANTU_NDEBUG auto old_size = send.size(); +#endif filter_list(keep_list, send); AKANTU_DEBUG_INFO("I had " << old_size << " elements to send to proc " << proc << " and " << keep_list.size() << " elements to keep. I have " << send.size() << " elements left."); } this->communicator.waitAll(send_requests); this->communicator.freeCommunicationRequest(send_requests); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::renumberElements( const ElementTypeMapArray & new_numbering) { for (auto && sr : iterate_send_recv) { for (auto && scheme_pair : communications.iterateSchemes(sr)) { auto & list = scheme_pair.second; for (auto && el : list) { if (new_numbering.exists(el.type, el.ghost_type)) el.element = new_numbering(el); } } } } /* -------------------------------------------------------------------------- */ UInt ElementSynchronizer::sanityCheckDataSize( const Array & elements, const SynchronizationTag &) const { UInt size = 0; size += sizeof(SynchronizationTag); // tag size += sizeof(UInt); // comm_desc.getNbData(); size += sizeof(UInt); // comm_desc.getProc(); size += sizeof(UInt); // mesh.getCommunicator().whoAmI(); // barycenters size += (elements.size() * mesh.getSpatialDimension() * sizeof(Real)); return size; } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::packSanityCheckData( CommunicationDescriptor & comm_desc) const { auto & buffer = comm_desc.getBuffer(); buffer << comm_desc.getTag(); buffer << comm_desc.getNbData(); buffer << comm_desc.getProc(); buffer << mesh.getCommunicator().whoAmI(); auto & send_element = comm_desc.getScheme(); /// pack barycenters in debug mode for (auto && element : send_element) { Vector barycenter(mesh.getSpatialDimension()); mesh.getBarycenter(element, barycenter); buffer << barycenter; } } /* -------------------------------------------------------------------------- */ void ElementSynchronizer::unpackSanityCheckData( CommunicationDescriptor & comm_desc) const { auto & buffer = comm_desc.getBuffer(); const auto & tag = comm_desc.getTag(); auto nb_data = comm_desc.getNbData(); auto proc = comm_desc.getProc(); auto rank = mesh.getCommunicator().whoAmI(); decltype(nb_data) recv_nb_data; decltype(proc) recv_proc; decltype(rank) recv_rank; SynchronizationTag t; buffer >> t; buffer >> recv_nb_data; buffer >> recv_proc; buffer >> recv_rank; AKANTU_DEBUG_ASSERT( t == tag, "The tag received does not correspond to the tag expected"); AKANTU_DEBUG_ASSERT( nb_data == recv_nb_data, "The nb_data received does not correspond to the nb_data expected"); AKANTU_DEBUG_ASSERT(UInt(recv_rank) == proc, "The rank received does not correspond to the proc"); AKANTU_DEBUG_ASSERT(recv_proc == UInt(rank), "The proc received does not correspond to the rank"); auto & recv_element = comm_desc.getScheme(); auto spatial_dimension = mesh.getSpatialDimension(); for (auto && element : recv_element) { Vector barycenter_loc(spatial_dimension); mesh.getBarycenter(element, barycenter_loc); Vector barycenter(spatial_dimension); buffer >> barycenter; auto dist = barycenter_loc.distance(barycenter); if (not Math::are_float_equal(dist, 0.)) AKANTU_EXCEPTION("Unpacking an unknown value for the element " << element << "(barycenter " << barycenter_loc << " != buffer " << barycenter << ") [" << dist << "] - tag: " << tag << " comm from " << proc << " to " << rank); } } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/test/test_model/patch_tests/CMakeLists.txt b/test/test_model/patch_tests/CMakeLists.txt index 81e811d25..d845ec350 100644 --- a/test/test_model/patch_tests/CMakeLists.txt +++ b/test/test_model/patch_tests/CMakeLists.txt @@ -1,64 +1,117 @@ #=============================================================================== # @file CMakeLists.txt # # @author David Simon Kammer # @author Nicolas Richart # # @date creation: Fri Oct 22 2010 # @date last modification: Mon Dec 07 2015 # # @brief configuration for patch tests # # @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 # #=============================================================================== add_subdirectory(data) register_gtest_sources(SOURCES patch_test_linear_elastic_explicit.cc PACKAGE solid_mechanics FILES_TO_COPY data/material_check_stress_plane_strain.dat data/material_check_stress_plane_stress.dat) register_gtest_sources(SOURCES patch_test_linear_elastic_implicit.cc PACKAGE solid_mechanics implicit FILES_TO_COPY data/material_check_stress_plane_strain.dat data/material_check_stress_plane_stress.dat) + register_gtest_sources(SOURCES patch_test_linear_anisotropic_explicit.cc - PACKAGE solid_mechanics - FILES_TO_COPY data/material_anisotropic.dat - UNSTABLE) + PACKAGE solid_mechanics lapack + FILES_TO_COPY data/material_anisotropic.dat) register_gtest_sources(SOURCES test_lumped_mass.cc PACKAGE solid_mechanics FILES_TO_COPY data/material_lumped.dat) register_gtest_sources( SOURCES patch_test_linear_heat_transfer_explicit.cc FILES_TO_COPY data/heat_transfer_input.dat PACKAGE heat_transfer) register_gtest_sources( SOURCES patch_test_linear_heat_transfer_static.cc FILES_TO_COPY data/heat_transfer_input.dat PACKAGE heat_transfer implicit) +register_gtest_sources( + SOURCES patch_test_linear_heat_transfer_implicit.cc + FILES_TO_COPY data/heat_transfer_input.dat + PACKAGE heat_transfer implicit) + register_gtest_test(patch_test_linear FILES_TO_COPY ${PATCH_TEST_MESHES}) + +register_test(test_linear_elastic_explicit_python + SCRIPT patch_test_linear_elastic_explicit.py + PYTHON + PACKAGE python_interface + FILES_TO_COPY patch_test_linear_fixture.py + FILES_TO_COPY patch_test_linear_solid_mechanics_fixture.py + FILES_TO_COPY ${PATCH_TEST_MESHES}) + +register_test(test_linear_elastic_static_python + SCRIPT patch_test_linear_elastic_static.py + PYTHON + PACKAGE python_interface implicit + FILES_TO_COPY patch_test_linear_fixture.py + FILES_TO_COPY patch_test_linear_solid_mechanics_fixture.py) + +register_test(test_linear_anisotropic_explicit_python + SCRIPT patch_test_linear_anisotropic_explicit.py + PYTHON + PACKAGE python_interface lapack + FILES_TO_COPY patch_test_linear_fixture.py + FILES_TO_COPY patch_test_linear_solid_mechanics_fixture.py + FILES_TO_COPY data/material_anisotropic.dat) + + +register_test(patch_test_linear_heat_transfer_explicit_python + SCRIPT patch_test_linear_heat_transfer_explicit.py + PYTHON + PACKAGE python_interface heat_transfer + FILES_TO_COPY patch_test_linear_fixture.py + FILES_TO_COPY patch_test_linear_heat_transfer_fixture.py + FILES_TO_COPY data/heat_transfer_input.dat) + +register_test(patch_test_linear_heat_transfer_static_python + SCRIPT patch_test_linear_heat_transfer_static.py + PYTHON + PACKAGE python_interface heat_transfer implicit + FILES_TO_COPY patch_test_linear_fixture.py + FILES_TO_COPY patch_test_linear_heat_transfer_fixture.py + FILES_TO_COPY data/heat_transfer_input.dat) + +register_test(patch_test_linear_heat_transfer_implicit_python + SCRIPT patch_test_linear_heat_transfer_implicit.py + PYTHON + PACKAGE python_interface heat_transfer implicit + FILES_TO_COPY patch_test_linear_fixture.py + FILES_TO_COPY patch_test_linear_heat_transfer_fixture.py + FILES_TO_COPY data/heat_transfer_input.dat) diff --git a/test/test_model/patch_tests/patch_test_linear_anisotropic_explicit.py b/test/test_model/patch_tests/patch_test_linear_anisotropic_explicit.py new file mode 100644 index 000000000..192ad4cf9 --- /dev/null +++ b/test/test_model/patch_tests/patch_test_linear_anisotropic_explicit.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python3 + +from patch_test_linear_solid_mechanics_fixture import TestPatchTestSMMLinear +import akantu +import unittest +import numpy as np + +# Stiffness tensor, rotated by hand +C = np.array([[[[112.93753505, 1.85842452538e-10, -4.47654358027e-10], + [1.85847317471e-10, 54.2334345331, -3.69840984824], + [-4.4764768395e-10, -3.69840984824, 56.848605217]], + [[1.85847781609e-10, 25.429294233, -3.69840984816], + [25.429294233, 3.31613847493e-10, -8.38797920011e-11], + [-3.69840984816, -8.38804581349e-11, -1.97875715813e-10]], + [[-4.47654358027e-10, -3.69840984816, 28.044464917], + [-3.69840984816, 2.09374961813e-10, 9.4857455224e-12], + [28.044464917, 9.48308098714e-12, -2.1367885239e-10]]], + [[[1.85847781609e-10, 25.429294233, -3.69840984816], + [25.429294233, 3.31613847493e-10, -8.38793479119e-11], + [-3.69840984816, -8.38795699565e-11, -1.97876381947e-10]], + [[54.2334345331, 3.31617400207e-10, 2.09372075233e-10], + [3.3161562385e-10, 115.552705733, -3.15093728886e-10], + [2.09372075233e-10, -3.15090176173e-10, 54.2334345333]], + [[-3.69840984824, -8.38795699565e-11, 9.48219280872e-12], + [-8.38795699565e-11, -3.1509195253e-10, 25.4292942335], + [9.48441325477e-12, 25.4292942335, 3.69840984851]]], + [[[-4.47653469848e-10, -3.69840984816, 28.044464917], + [-3.69840984816, 2.09374073634e-10, 9.48752187924e-12], + [28.044464917, 9.48552347779e-12, -2.1367885239e-10]], + [[-3.69840984824, -8.3884899027e-11, 9.48219280872e-12], + [-8.3884899027e-11, -3.150972816e-10, 25.4292942335], + [9.48041645188e-12, 25.4292942335, 3.69840984851]], + [[56.848605217, -1.97875493768e-10, -2.13681516925e-10], + [-1.97877270125e-10, 54.2334345333, 3.69840984851], + [-2.13683293282e-10, 3.69840984851, 112.93753505]]]]) + + +def foo(self): + + self.initModel( + akantu.SolidMechanicsModelOptions(akantu._explicit_lumped_mass), + "material_anisotropic.dat") + + coordinates = self.mesh.getNodes() + displacement = self.model.getDisplacement() + + # set the position of all nodes to the static solution + self.setLinearDOF(displacement, coordinates) + + for s in range(0, 100): + self.model.solveStep() + + ekin = self.model.getEnergy("kinetic") + self.assertAlmostEqual(0, ekin, delta=1e-16) + + self.checkDisplacements() + self.checkStrains() + + def foo(pstrain): + strain = (pstrain + pstrain.transpose()) / 2. + stress = np.zeros((self.dim, self.dim)) + + for i in range(0, self.dim): + for j in range(0, self.dim): + stress[i, j] = 0 + for k in range(0, self.dim): + for l in range(0, self.dim): + stress[i, j] += C[i][j][k][l] * strain(k, l) + return stress + self.checkStresses(foo) + + +suite = TestPatchTestSMMLinear.TYPED_TEST(foo, "AnisotropicExplicit") +unittest.TextTestRunner(verbosity=1).run(suite) diff --git a/test/test_model/patch_tests/patch_test_linear_elastic_explicit.py b/test/test_model/patch_tests/patch_test_linear_elastic_explicit.py new file mode 100644 index 000000000..979c3dc0f --- /dev/null +++ b/test/test_model/patch_tests/patch_test_linear_elastic_explicit.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python3 + +from patch_test_linear_solid_mechanics_fixture import TestPatchTestSMMLinear +import akantu + + +def foo(self): + + filename = "material_check_stress_plane_stress.dat" + if self.plane_strain: + filename = "material_check_stress_plane_strain.dat" + + self.initModel( + akantu.SolidMechanicsModelOptions(akantu._explicit_lumped_mass), + filename) + + coordinates = self.mesh.getNodes() + displacement = self.model.getDisplacement() + # set the position of all nodes to the static solution + self.setLinearDOF(displacement, coordinates) + + for s in range(0, 100): + self.model.solveStep() + + ekin = self.model.getEnergy("kinetic") + self.assertAlmostEqual(0, ekin, 1e-16) + self.checkAll() + + +TestPatchTestSMMLinear.TYPED_TEST(foo, "Explicit") diff --git a/test/test_model/patch_tests/patch_test_linear_elastic_static.py b/test/test_model/patch_tests/patch_test_linear_elastic_static.py new file mode 100644 index 000000000..d589dac21 --- /dev/null +++ b/test/test_model/patch_tests/patch_test_linear_elastic_static.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python3 + +from patch_test_linear_solid_mechanics_fixture import TestPatchTestSMMLinear +import akantu + + +def foo(self): + + filename = "material_check_stress_plane_stress.dat" + if self.plane_strain: + filename = "material_check_stress_plane_strain.dat" + + self.initModel(akantu.SolidMechanicsModelOptions(akantu._static), filename) + + solver = self.model.getNonLinearSolver() + solver.set("max_iterations", 2) + solver.set("threshold", 2e-4) + solver.set("convergence_type", akantu._scc_residual) + + self.model.solveStep() + + self.checkAll() + + +TestPatchTestSMMLinear.TYPED_TEST(foo, "Static") diff --git a/test/test_model/patch_tests/patch_test_linear_fixture.py b/test/test_model/patch_tests/patch_test_linear_fixture.py new file mode 100644 index 000000000..765790449 --- /dev/null +++ b/test/test_model/patch_tests/patch_test_linear_fixture.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python3 + +import akantu +import unittest +import numpy as np +import sys + + +class TestPatchTestLinear(unittest.TestCase): + + alpha = np.array([[0.01, 0.02, 0.03, 0.04], + [0.05, 0.06, 0.07, 0.08], + [0.09, 0.10, 0.11, 0.12]]) + + gradient_tolerance = 1e-13 + result_tolerance = 1e-14 + dofs_tolerance = 1e-15 + + def __init__(self, test_name, elem_type_str, functor=None): + + self.test_name = test_name + self.functor = functor + self.elem_type = akantu.getElementTypes()[elem_type_str] + self.elem_type_str = elem_type_str + self.dim = akantu.Mesh.getSpatialDimension(self.elem_type) + + super().__init__(test_name) + + def _internal_call(self): + self.functor(self) + + def __getattr__(self, key): + if key == self.test_name: + return self._internal_call + + def setUp(self): + self.mesh = akantu.Mesh(self.dim, self.elem_type_str) + self.mesh.read(str(self.elem_type_str) + ".msh") + akantu.MeshUtils.buildFacets(self.mesh) + self.mesh.createBoundaryGroupFromGeometry() + self.model = self.model_type(self.mesh) + + def tearDown(self): + del self.model + del self.mesh + + def initModel(self, method, material_file): + akantu.initialize(material_file) + akantu.setDebugLevel(akantu.dblError) + self.model.initFull(method) + self.applyBC() + + def applyBC(self): + boundary = self.model.getBlockedDOFs() + for name, eg in self.mesh.getElementGroups().items(): + nodes = eg['node_group'] + boundary[nodes, :] = True + + def applyBConDOFs(self, dofs): + coordinates = self.mesh.getNodes() + + for name, eg in self.mesh.getElementGroups().items(): + nodes = eg['node_group'].flatten() + dofs[nodes] = self.setLinearDOF(dofs[nodes], + coordinates[nodes]) + + def prescribed_gradient(self, dof): + gradient = self.alpha[:dof.shape[1], 1:self.dim + 1] + return gradient + + def checkGradient(self, gradient, dofs): + pgrad = self.prescribed_gradient(dofs).T + gradient = gradient.reshape(gradient.shape[0], *pgrad.shape) + diff = gradient[:] - pgrad + norm = np.abs(pgrad).max() + gradient_error = np.abs(diff).max() / norm + self.assertAlmostEqual(0, gradient_error, + delta=self.gradient_tolerance) + + def checkResults(self, presult_func, results, dofs): + presult = presult_func(self.prescribed_gradient(dofs)).flatten() + remaining_size = np.prod(np.array(results.shape[1:])) + results = results.reshape((results.shape[0], remaining_size)) + + for result in results: + diff = result - presult + norm = np.abs(result).max() + if norm == 0: + result_error = np.abs(diff).max() + else: + result_error = np.abs(diff).max() / norm + + self.assertAlmostEqual(0., result_error, + delta=self.result_tolerance) + + def setLinearDOF(self, dof, coord): + nb_dofs = dof.shape[1] + dof[:] = np.einsum('ik,ak->ai', + self.alpha[:nb_dofs, 1:self.dim + 1], coord) + for i in range(0, nb_dofs): + dof[:, i] += self.alpha[i, 0] + + return dof + + def checkDOFs(self, dofs): + coordinates = self.mesh.getNodes() + ref_dofs = np.zeros_like(dofs) + self.setLinearDOF(ref_dofs, coordinates) + diff = dofs - ref_dofs + dofs_error = np.abs(diff).max() + self.assertAlmostEqual(0., dofs_error, delta=self.dofs_tolerance) + + @classmethod + def TYPED_TEST(cls, functor, label): + suite = unittest.TestSuite() + + for type_name, _type in akantu.getElementTypes().items(): + if type_name == "_point_1": + continue + + method_name = type_name + '_' + label + test_case = cls(method_name, type_name, functor) + suite.addTest(test_case) + + result = unittest.TextTestRunner(verbosity=1).run(suite) + ret = not result.wasSuccessful() + sys.exit(ret) diff --git a/test/test_model/patch_tests/patch_test_linear_heat_transfer_explicit.py b/test/test_model/patch_tests/patch_test_linear_heat_transfer_explicit.py new file mode 100644 index 000000000..e0f89c671 --- /dev/null +++ b/test/test_model/patch_tests/patch_test_linear_heat_transfer_explicit.py @@ -0,0 +1,24 @@ +#!/usr/bin/env python3 + +from patch_test_linear_heat_transfer_fixture import TestPatchTestHTMLinear +import akantu + + +def foo(self): + + self.initModel( + akantu.HeatTransferModelOptions(akantu._explicit_lumped_mass), + "heat_transfer_input.dat") + + coordinates = self.mesh.getNodes() + temperature = self.model.getTemperature() + # set the position of all nodes to the static solution + self.setLinearDOF(temperature, coordinates) + + for s in range(0, 100): + self.model.solveStep() + + self.checkAll() + + +TestPatchTestHTMLinear.TYPED_TEST(foo, "Explicit") diff --git a/test/test_model/patch_tests/patch_test_linear_heat_transfer_fixture.py b/test/test_model/patch_tests/patch_test_linear_heat_transfer_fixture.py new file mode 100644 index 000000000..de08c4b0d --- /dev/null +++ b/test/test_model/patch_tests/patch_test_linear_heat_transfer_fixture.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python3 + +import patch_test_linear_fixture +import akantu + + +class TestPatchTestHTMLinear(patch_test_linear_fixture.TestPatchTestLinear): + + model_type = akantu.HeatTransferModel + + def applyBC(self): + super().applyBC() + temperature = self.model.getTemperature() + self.applyBConDOFs(temperature) + + def checkAll(self): + temperature = self.model.getTemperature() + C = self.model.getParamMatrix("conductivity") + self.checkDOFs(temperature) + self.checkGradient(self.model.getTemperatureGradient(self.elem_type), + temperature) + + self.prescribed_gradient(temperature) + self.checkResults(lambda grad_T: C.dot(grad_T.T), + self.model.getKgradT(self.elem_type), + temperature) diff --git a/test/test_model/patch_tests/patch_test_linear_heat_transfer_implicit.cc b/test/test_model/patch_tests/patch_test_linear_heat_transfer_implicit.cc index c953048cf..c14e0207b 100644 --- a/test/test_model/patch_tests/patch_test_linear_heat_transfer_implicit.cc +++ b/test/test_model/patch_tests/patch_test_linear_heat_transfer_implicit.cc @@ -1,24 +1,21 @@ /* -------------------------------------------------------------------------- */ #include "patch_test_linear_heat_transfer_fixture.hh" /* -------------------------------------------------------------------------- */ -TYPED_TEST(TestPatchTestHTMLinear, Explicit) { - this->initModel(_explicit_lumped_mass, "heat_transfer_input.dat"); +TYPED_TEST(TestPatchTestHTMLinear, Implicit) { + this->initModel(_implicit_dynamic, "heat_transfer_input.dat"); const auto & coordinates = this->mesh->getNodes(); auto & temperature = this->model->getTemperature(); // set the position of all nodes to the static solution for (auto && tuple : zip(make_view(coordinates, this->dim), make_view(temperature, 1))) { this->setLinearDOF(std::get<1>(tuple), std::get<0>(tuple)); } for (UInt s = 0; s < 100; ++s) { this->model->solveStep(); } - auto eth = this->model->getEnergy("thermal"); - EXPECT_NEAR(0, eth, 1e-16); - this->checkAll(); } diff --git a/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.py b/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.py new file mode 100644 index 000000000..25c4e61af --- /dev/null +++ b/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.py @@ -0,0 +1,151 @@ +#!/usr/bin/env python3 + +import patch_test_linear_fixture +import numpy as np +import akantu + + +# custom material (this patch test also checks for custom material features) +class LocalElastic: + + # declares all the internals + def initMaterial(self, internals, params): + self.E = params['E'] + self.nu = params['nu'] + self.rho = params['rho'] + # First Lame coefficient + self.lame_lambda = self.nu * self.E / ( + (1. + self.nu) * (1. - 2. * self.nu)) + # Second Lame coefficient (shear modulus) + self.lame_mu = self.E / (2. * (1. + self.nu)) + + # declares all the internals + @staticmethod + def registerInternals(): + return ['potential'] + + # declares all the internals + @staticmethod + def registerInternalSizes(): + return [1, 1] + + # declares all the parameters that could be parsed + @staticmethod + def registerParam(): + return ['E', 'nu'] + + # declares all the parameters that are needed + def getPushWaveSpeed(self, params): + return np.sqrt((self.lame_lambda + 2 * self.lame_mu) / self.rho) + + # compute small deformation tensor + @staticmethod + def computeEpsilon(grad_u): + return 0.5 * (grad_u + np.einsum('aij->aji', grad_u)) + + # constitutive law + def computeStress(self, grad_u, sigma, internals, params): + n_quads = grad_u.shape[0] + grad_u = grad_u.reshape((n_quads, self.dim, self.dim)) + epsilon = self.computeEpsilon(grad_u) + sigma = sigma.reshape((n_quads, self.dim, self.dim)) + + trace = np.einsum('aii->a', grad_u) + + if self.dim == 1: + sigma[:, :, :] = self.E * epsilon + + else: + sigma[:, :, :] = ( + np.einsum('a,ij->aij', trace, + self.lame_lambda * np.eye(self.dim)) + + 2.*self.lame_mu * epsilon) + + + # constitutive law tangent modulii + def computeTangentModuli(self, grad_u, tangent, internals, params): + + n_quads = tangent.shape[0] + tangent = tangent.reshape(n_quads, 3, 3) + + Miiii = self.lame_lambda + 2 * self.lame_mu + Miijj = self.lame_lambda + Mijij = self.lame_mu + + tangent[:, 0, 0] = Miiii + tangent[:, 1, 1] = Miiii + tangent[:, 0, 1] = Miijj + tangent[:, 1, 0] = Miijj + tangent[:, 2, 2] = Mijij + + # computes the energy density + def getEnergyDensity(self, energy_type, energy_density, + grad_u, stress, internals, params): + + nquads = stress.shape[0] + stress = stress.reshape(nquads, 2, 2) + grad_u = grad_u.reshape((nquads, 2, 2)) + + if energy_type != 'potential': + raise RuntimeError('not known energy') + + epsilon = self.computeEpsilon(grad_u) + + energy_density[:, 0] = ( + 0.5 * np.einsum('aij,aij->a', stress, epsilon)) + + +class TestPatchTestSMMLinear(patch_test_linear_fixture.TestPatchTestLinear): + + plane_strain = True + model_type = akantu.SolidMechanicsModel + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + def initModel(self, method, material_file): + mat.__dict__['dim'] = self.dim + super().initModel(method, material_file) + + def applyBC(self): + super().applyBC() + displacement = self.model.getDisplacement() + self.applyBConDOFs(displacement) + + def checkAll(self): + displacement = self.model.getDisplacement() + mat = self.model.getMaterial(0) + + self.checkDOFs(displacement) + self.checkGradient(mat.getGradU(self.elem_type), displacement) + + def foo(pstrain): + nu = self.model.getMaterial(0).getParamReal("nu") + E = self.model.getMaterial(0).getParamReal("E") + + strain = (pstrain + pstrain.transpose()) / 2. + trace = strain.trace() + + _lambda = nu * E / ((1 + nu) * (1 - 2 * nu)) + mu = E / (2 * (1 + nu)) + + if (not self.plane_strain): + _lambda = nu * E / (1 - nu * nu) + + stress = np.zeros((self.dim, self.dim)) + + if self.dim == 1: + stress[0, 0] = E * strain[0, 0] + else: + stress[:, :] = ( + _lambda * trace * np.eye(self.dim) + 2 * mu * strain) + return stress + + self.checkResults(foo, + mat.getStress(self.elem_type), + displacement) + + + +mat = LocalElastic() +akantu.registerNewPythonMaterial(mat, "local_elastic")