Diffusion akantu (bugfixes/contact_mechanics_model) Browse

contact_mechan…

custom_material.py
View Options

			#!/usr/bin/env python3
			# pylint: disable=missing-module-docstring
			# pylint: disable=missing-function-docstring

			import numpy as np
			import akantu as aka

			spatial_dimension = 2

			# ------------------------------------------------------------------------------
			class LinearCohesive(aka.MaterialCohesive):
			"""Material Cohesive Linear."""

			def __init__(self, model, _id):
			super().__init__(model, _id)
			super().registerParamReal(
			"G_c", aka._pat_readable \| aka._pat_parsable, "Fracture energy"
			)
			super().registerParamReal("beta", aka._pat_readable \| aka._pat_parsable, "beta")
			self.registerInternalReal("delta_max", 1)
			self.beta = 0
			self.sigma_c = 0
			self.delta_c = 0

			def initMaterial(self):
			super().initMaterial()
			self.sigma_c = self.getReal("sigma_c")
			self.G_c = self.getReal("G_c")
			self.beta = self.getReal("beta")
			self.delta_c = 2 * self.G_c / self.sigma_c

			def checkInsertion(self, check_only):
			model = self.getModel()
			facets = self.getFacetFilter()
			inserter = model.getElementInserter()

			for type_facet in facets.elementTypes(dim=(spatial_dimension - 1)):
			facet_filter = facets(type_facet)
			nb_facet = facet_filter.shape[0]
			if nb_facet == 0:
			continue

			fe_engine = model.getFEEngine("FacetsFEEngine")

			facets_check = inserter.getCheckFacets(type_facet)
			insertion = inserter.getInsertionFacets(type_facet)

			nb_quad_facet = fe_engine.getNbIntegrationPoints(type_facet)
			normals = fe_engine.getNormalsOnIntegrationPoints(type_facet)
			facets_stresses = model.getStressOnFacets(type_facet).reshape(
			normals.shape[0] // nb_quad_facet,
			nb_quad_facet,
			2,
			spatial_dimension,
			spatial_dimension,
			)

			tangents = model.getTangents(type_facet)

			for facet, facet_stresses in zip(facet_filter, facets_stresses):
			if not facets_check[facet]:
			continue

			ref_stress = 0

			for q, quad_stresses in enumerate(facet_stresses):
			current_quad = facet * nb_quad_facet + q
			normal = normals[current_quad, :].ravel()
			tangent = tangents[current_quad, :].ravel()

			stress_1 = quad_stresses.T[0]
			stress_2 = quad_stresses.T[1]

			avg_stress = stress_1 + stress_2 / 2.0
			traction = avg_stress.dot(normal)

			n = traction.dot(normal)
			n = max(0, n)
			t = traction.dot(tangent)

			ref_stress = max(
			ref_stress, np.sqrt(n * n + t * t / (self.beta**2))
			)

			if ref_stress > self.sigma_c:
			insertion[facet] = True

			# constitutive law
			def computeTraction(self, normals, el_type, ghost_type):
			openings = self.getOpening(el_type, ghost_type)
			tractions = self.getTraction(el_type, ghost_type)

			delta_max = self.getInternalReal("delta_max")(el_type)

			for el in range(normals.shape[0]):
			normal = normals[el].ravel()
			opening = openings[el].ravel()

			delta_n = opening.dot(normal) * normal
			delta_s = opening - delta_n

			delta = (
			self.beta * np.linalg.norm(delta_s) 2 + np.linalg.norm(delta_n) 2
			)

			delta_max[el] = max(delta_max[el], delta)

			tractions[el, :] = (
			(delta * delta_s + delta_n)
			* self.sigma_c
			/ delta
			* (1 - delta / self.delta_c)
			)


			# register material to the MaterialFactory
			def allocator(_dim, unused, model, _id):
			return LinearCohesive(model, _id)


			mat_factory = aka.MaterialFactory.getInstance()
			mat_factory.registerAllocator("local_cohesive", allocator)

			# -------------------------------------------------------------------------
			# Initialization
			# -------------------------------------------------------------------------
			aka.parseInput("local_material.dat")
			mesh = aka.Mesh(spatial_dimension)
			mesh.read("plate.msh")

			model = aka.SolidMechanicsModelCohesive(mesh)
			model.initFull(_analysis_method=aka._static, _is_extrinsic=True)
			model.initNewSolver(aka._explicit_lumped_mass)

			model.setBaseName("plate")
			model.addDumpFieldVector("displacement")
			model.addDumpFieldVector("external_force")
			model.addDumpField("strain")
			model.addDumpField("stress")
			model.addDumpField("blocked_dofs")

			model.setBaseNameToDumper("cohesive elements", "cohesive")
			model.addDumpFieldVectorToDumper("cohesive elements", "displacement")
			model.addDumpFieldToDumper("cohesive elements", "damage")
			model.addDumpFieldVectorToDumper("cohesive elements", "tractions")
			model.addDumpFieldVectorToDumper("cohesive elements", "opening")

			# -------------------------------------------------------------------------
			# Boundary conditions
			# -------------------------------------------------------------------------
			model.applyBC(aka.FixedValue(0.0, aka._x), "XBlocked")
			model.applyBC(aka.FixedValue(0.0, aka._y), "YBlocked")

			trac = np.zeros(spatial_dimension)
			traction = 0.095
			trac[int(aka._y)] = traction

			model.getExternalForce()[:] = 0
			model.applyBC(aka.FromTraction(trac), "Traction")

			print("Solve for traction ", traction)
			solver = model.getNonLinearSolver("static")
			solver.set("max_iterations", 100)
			solver.set("threshold", 1e-10)
			solver.set("convergence_type", aka.SolveConvergenceCriteria.residual)

			model.solveStep("static")
			model.dump()
			model.dump("cohesive elements")

			model.setTimeStep(model.getStableTimeStep() * 0.1)

			maxsteps = 100

			for i in range(0, maxsteps):
			print("{0}/{1}".format(i, maxsteps))
			model.checkCohesiveStress()
			model.solveStep("explicit_lumped")
			if i % 10 == 0:
			model.dump()
			model.dump("cohesive elements")

rAKA/examples/python/cohesive/custom_material.py3b7d23e08456bugfixes/contact_mechan…

custom_material.pyView Options

rAKA/examples/python/cohesive/custom_material.py
3b7d23e08456bugfixes/contact_mechan…

custom_material.py
View Options