# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
import numpy as np
from scipy.sparse import coo_matrix, diags
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.MeshContainers.Filters.FilterOperators import UnionFilter, IntersectionFilter, FilterLike
from Muscat.MeshTools.MeshCreationTools import CreateMeshOf
from Muscat.FE.IntegrationRules import LagrangeIsoParamQuadrature
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
from Muscat.LinAlg.Transform import Transform
from Muscat.FE.KR.KRBase import KRBaseVector
from Muscat.FE.Fields.FEField import FEField
from Muscat.MeshContainers.Mesh import Mesh
from Muscat.Helpers.Decorators import froze_it
from Muscat.Helpers.ProgressBar import PrintProgressBar
[docs]
@froze_it
class KRMortar(KRBaseVector):
def __init__(self):
super().__init__()
self.type = "Mortar"
self.useSurface = "first_surface" # [ "mean_surface", "first_surface", "second_surface","flat"]
self.onII =[]
self.ang = 30/180.*np.pi # tolerance angle
self.originSystem.keepOrthogonal = True
self.originSystem.keepNormalized = True
self.targetSystem.keepOrthogonal = True
self.targetSystem.keepNormalized = True
self._debug_IntegrationMesh = None
self.weights_IPoints = None
self.meshIOps = None
self.meshIIOps = None
self.meshI_IPoints = None
self.meshII_IPoints = None
[docs]
def From(self,offset=None,first=None,second=None):
if offset is not None:
self.originSystem.SetOffset(offset)
if first is not None:
self.originSystem.SetFirst(first)
if second is not None:
self.originSystem.SetSecond(second)
return self
[docs]
def To(self,offset=None,first=None,second=None):
if offset is not None:
self.targetSystem.SetOffset(offset)
if first is not None:
self.targetSystem.SetFirst(first)
if second is not None:
self.targetSystem.SetSecond(second)
return self
[docs]
def SideI(self,zone):
return self.On(zone)
[docs]
def SideII(self,zone):
if type(zone) is list:
self.onII.extend(zone)
else:
self.onII.append(zone)
self.onII = list(set(self.onII))
return self
def _GetOnIIElementFilter(self) -> FilterLike:
"""Create a filter with the data from onII
Returns
-------
FilterLike
the filter to select the elements to treat
"""
eTags = []
filters = []
for obj in self.onII:
if isinstance(obj, str):
eTags.append(obj)
else:
filters.append(obj)
if len(eTags):
filters.append(ElementFilter(eTag=eTags))
return UnionFilter(filters=filters)
[docs]
def GenerateEquations(self,meshI,fields,CH=None,meshII=None,fieldsII=None):
CH = self._GetConstraintHolder(CH)
if meshII is None:
meshII = meshI
fieldsII = fields
fieldDicI = {f.name:f for f in fields }
fieldDicII = {f.name:f for f in fields }
offsetsI , __fieldOffsetsI, totalNumberOfDofsI = self._ComputeOffsets(fields)
offsetsII, __fieldOffsetsII, totalNumberOfDofsII = self._ComputeOffsets(fieldsII)
def AddElementToViz(visualizationMesh, points, t,tag=None): # pragma: no cover
if visualizationMesh is None:
return
n = visualizationMesh.GetNumberOfNodes()
visualizationMesh.SetNodes(np.concatenate((visualizationMesh.nodes, points), axis=0),generateOriginalIDs=True)
data = visualizationMesh.GetElementsOfType(t)
cpt = data.AddNewElement(range(n,n+data.GetNumberOfNodesPerElement()),0)
if tag is not None:
data.GetTag(tag).AddToTag(cpt-1)
bar_p = LagrangeIsoParamQuadrature[ED.Bar_2].points
bar_w = LagrangeIsoParamQuadrature[ED.Bar_2].weights
barSpaceIPValues = LagrangeSpaceGeo[ED.Bar_2].GetSpaceEvaluatedAt(bar_p,bar_w)
tri_p = LagrangeIsoParamQuadrature[ED.Triangle_3].points
tri_w = LagrangeIsoParamQuadrature[ED.Triangle_3].weights
triSpaceIPValues = LagrangeSpaceGeo[ED.Triangle_3].GetSpaceEvaluatedAt(tri_p,tri_w)
tet_p = LagrangeIsoParamQuadrature[ED.Tetrahedron_4].points
tet_w = LagrangeIsoParamQuadrature[ED.Tetrahedron_4].weights
tetSpaceIPValues = LagrangeSpaceGeo[ED.Tetrahedron_4].GetSpaceEvaluatedAt(tet_p,tet_w)
meshI_IPoints = []
meshII_IPoints = []
weights_IPoints = []
mesh_IElement = []
mesh_IIElement = []
# reference systems
OS = self.originSystem.GetOrthoNormalBase()
TS = self.targetSystem.GetOrthoNormalBase()
# working memory
# for bars
integrationBar = np.zeros((2,3))
iwBar = np.zeros_like(bar_w)
ipBar = np.zeros((len(bar_w),3))
# for triangles
integrationTri = np.zeros((3,3))
iwTri = np.zeros_like(tri_w)
ipTri = np.zeros((len(tri_w),3))
# for tetra
integrationTet = np.zeros((4,3))
iwTet = np.zeros_like(tet_w)
ipTet = np.zeros((len(tet_w),3))
# loop over every element in the two meshes
#computation of the integration points
for selection1 in self._GetOnElementFilter()(meshI):
cpt = 0
for elementId1 in selection1.indices:
PrintProgressBar(cpt,len(selection1.indices))
cpt += 1
#pos of node in the real spaces
_nodes1 = meshI.nodes[selection1.elements.connectivity[elementId1],:]
#normal in the real space
if ED.dimensionality[selection1.elementType] < 3:
_normal1,_barycenter1, inCellVector1 = self.__ComputeNormal(meshI,selection1.elementType,selection1.elements,elementId1)
#barycenter in the projection space
#barycenter1 = OS.ApplyTransform(_barycenter1)
#normal in the projection space
normal1 = OS.ApplyTransformDirection(_normal1)
#nodes in the projection space
nodes1 = OS.ApplyTransform(_nodes1)
#bounding box
min1 = np.min(nodes1,axis=0)
max1 = np.max(nodes1,axis=0)
for selection2 in IntersectionFilter((ElementFilter(dimensionality=ED.dimensionality[selection1.elementType]),self._GetOnIIElementFilter())) (meshII):
for elementId2 in selection2.indices:
_nodes2 = meshII.nodes[selection2.elements.connectivity[elementId2],:]
nodes2 = TS.ApplyTransform(_nodes2)
# check if bounding box intersect
min2 = np.min(nodes2,axis=0)
max2 = np.max(nodes2,axis=0)
intersection = True
tol = 0.1*max( np.max(max1-min1),np.max(max2-min2))
for i in range(len(min1)):
intersection = intersection and ((min2[i] <= max1[i]+tol and min2[i] >= min1[i]-tol) or (
max2[i] <= max1[i]+tol and max2[i] >= min1[i]-tol) or (
min1[i] <= max2[i]+tol and min1[i] >= min2[i]-tol) or (
max1[i] <= max2[i]+tol and max1[i] >= min2[i]-tol))
if not intersection :
continue
if ED.dimensionality[selection2.elementType] < 3:
_normal2,_barycenter2, inCellVector2 = self.__ComputeNormal(meshII,selection2.elementType,selection2.elements,elementId2)
normal2 = TS.ApplyTransformDirection(_normal2)
# only for 1d and 2d elements
# check normal are aligned using the self.ang
# the orientation is not important
angle = np.arccos(normal1.dot(normal2))
if abs(angle) > self.ang and abs(angle-np.pi) > self.ang:
continue
# we have intersection of 2 elements
if ED.dimensionality[selection2.elementType] == 1:
# element of dimension 1 (bars)
if self.useSurface == "mean_surface":
lineNormal = (normal1-normal2)/2
lineNormal /= np.linalg.norm(lineNormal)
base = np.array([-lineNormal[1], lineNormal[0], lineNormal[2]])
offset = (min1+min2+max1+max2)/4
# projection of the point in this line
elif self.useSurface == "first_surface":
base = np.array([-normal1[1],normal1[0],normal1[2]])
offset = (min1+max1)/2
elif self.useSurface == "second_surface":
base = np.array([-normal2[1],normal2[0],normal2[2]])
offset = (min2+max2)/2
else:
raise(Exception(f"Method {self.useSurface} not available!!!"))
second = base*1
for i,v in enumerate(second):
if v == 0:
second[i] += 1
break
else:
second[0] += 1
T = Transform(offset=offset,first=base, second=second )
proj1_ = T.ApplyTransform(nodes1)
proj2_ = T.ApplyTransform(nodes2)
proj1 = proj1_[:,0]
proj2 = proj2_[:,0]
#compute intersection
projMin = max(min(proj1),min(proj2))
projMax = min(max(proj1),max(proj2))
if projMax < projMin:
continue
integrationBar[0,0 ] = projMin
integrationBar[1,0 ] = projMax
# compute the coordinate of the integration points
for ip_nb in range(len(bar_w)):
Jacobian, JDet, JInv = barSpaceIPValues.GetJacobianDeterminantAndInverseAtIP(ip_nb,integrationBar)
iwBar[ip_nb] = JDet*bar_w[ip_nb]
ipBar[ip_nb,:] = np.dot(barSpaceIPValues.valN[ip_nb],integrationBar)
# if the integration bar is degenerated we skip it
if JDet < 1e-8:
continue
# transfer the coordinate to the original meshes (1 and 2)
integration_coords_II = T.ApplyInvTransform(ipBar)
## compute integration point in the meshI, meshII
original_int_points = OS.ApplyInvTransform(integration_coords_II)
target_int_points = TS.ApplyInvTransform(integration_coords_II)
if self._debug_IntegrationMesh is not None: # pragma: no cover
AddElementToViz(self._debug_IntegrationMesh,original_int_points,ED.Bar_2,"eInt_I")
AddElementToViz(self._debug_IntegrationMesh,target_int_points,ED.Bar_2,"eInt_II")
# Append points and weights
weights_IPoints.extend(iwBar)
mesh_IElement.extend(np.full(len(bar_w),elementId1))
mesh_IIElement.extend(np.full(len(bar_w),elementId2))
for ip_nb in range(len(bar_w)):
meshI_IPoints.append(original_int_points[ip_nb,:])
meshII_IPoints.append(target_int_points[ip_nb,:])
if self._debug_IntegrationMesh is not None:# pragma: no cover
AddElementToViz(self._debug_IntegrationMesh,original_int_points[ip_nb:ip_nb+1,:],ED.Point_1,"int_I")
AddElementToViz(self._debug_IntegrationMesh,target_int_points[ip_nb:ip_nb+1,:],ED.Point_1,"int_II")
elif ED.dimensionality[selection2.elementType] == 2:
# element of dimension 2 (triangles)
if self.useSurface == "flat":
lineNormal = [0,0,1]
elif self.useSurface == "mean_surface":
lineNormal = (normal1-normal2)/2
lineNormal /= np.linalg.norm(lineNormal)
# projection of the point in this line
elif self.useSurface == "first_surface":
lineNormal = normal1
elif self.useSurface == "second_surface":
lineNormal = normal2
#T = Transform(first=[-lineNormal[1], lineNormal[0], lineNormal[2]],second=[-lineNormal[2], lineNormal[1], lineNormal[0]])
T = Transform()
T.SetOpUsingThird(lineNormal)
proj1 = T.ApplyTransform(nodes1)
proj2 = T.ApplyTransform(nodes2)
inter = Intersection(proj1*[1,1,0],range(selection1.elements.GetNumberOfNodesPerElement()),
proj2*[1,1,0],range(selection2.elements.GetNumberOfNodesPerElement()),tol/100000.)
if len(inter) < 4 :
# insufficient point to build as 2D domain
continue
if len(inter) == 4 :
# we have a triangle (first point is repeated at the end)
center = inter[0,:]
inter = inter[1:3,:]
else:
center = np.sum(inter[0:-1,:],axis=0)/(inter.shape[0]-1)
# we use the center as the first point for all the triangles
integrationTri[0,:] = center
for i in range(inter.shape[0]):
# copy the coord to generate a triangle
integrationTri[1:3,:] = inter[i:i+2,:]
# compute the coordinate of the integration points
for ip_nb in range(len(tri_w)):
Jacobian, JDet, JInv = triSpaceIPValues.GetJacobianDeterminantAndInverseAtIP(ip_nb,integrationTri)
iwTri[ip_nb] = JDet*tri_w[ip_nb]
ipTri[ip_nb,:] = np.dot(triSpaceIPValues.valN[ip_nb],integrationTri)
# if the integration triangle is degenerated we skip it
if JDet < 1e-8:
continue
# transfer the coordinate to the original meshes (1 and 2)
integration_coords_II = T.ApplyInvTransform(ipTri)
if self._debug_IntegrationMesh is not None:# pragma: no cover
AddElementToViz(self._debug_IntegrationMesh,OS.ApplyInvTransform(T.ApplyInvTransform(integrationTri)),ED.Triangle_3,"eInt_I")
AddElementToViz(self._debug_IntegrationMesh,TS.ApplyInvTransform(T.ApplyInvTransform(integrationTri)),ED.Triangle_3,"eInt_II")
## compute integration point in the meshI, meshII
original_int_points = OS.ApplyInvTransform(integration_coords_II)
target_int_points = TS.ApplyInvTransform(integration_coords_II)
if self._debug_IntegrationMesh is not None:# pragma: no cover
AddElementToViz(self._debug_IntegrationMesh,original_int_points,ED.Triangle_3,"eInt_I")
AddElementToViz(self._debug_IntegrationMesh,target_int_points,ED.Triangle_3,"eInt_II")
# Append points and weights
weights_IPoints.extend(iwTri)
mesh_IElement.extend(np.full(len(tri_w),elementId1))
mesh_IIElement.extend(np.full(len(tri_w),elementId2))
for ip_nb in range(len(tri_w)):
meshI_IPoints.append(original_int_points[ip_nb,:])
meshII_IPoints.append(target_int_points[ip_nb,:])
AddElementToViz(self._debug_IntegrationMesh,original_int_points[ip_nb:ip_nb+1,:],ED.Point_1,"int_I")
AddElementToViz(self._debug_IntegrationMesh,target_int_points[ip_nb:ip_nb+1,:],ED.Point_1,"int_II")
elif ED.dimensionality[selection2.elementType] == 3:
status, points, tets = IntersectionOf2ConvexHulls(nodes1,nodes2)
if status == False:
continue
for i in range(tets.shape[0]):
# copy the coord to generate a tet
integrationTet[:,:] = points[tets[i,:],:]
for ip_nb in range(len(tet_w)):
Jacobian, JDet, JInv = tetSpaceIPValues.GetJacobianDeterminantAndInverseAtIP(ip_nb,integrationTet)
iwTet[ip_nb] = abs(JDet)*tet_w[ip_nb]
ipTet[ip_nb,:] = np.dot(tetSpaceIPValues.valN[ip_nb],integrationTet)
# if the integration triangle is degenerated we skip it
if abs(JDet) < 1e-8:
continue
# transfer the coordinate to the original meshes (1 and 2)
if self._debug_IntegrationMesh is not None: # pragma: no cover
AddElementToViz(self._debug_IntegrationMesh,OS.ApplyInvTransform(integrationTet),ED.Tetrahedron_4,"eInt_I_"+str(elementId1)+"_"+str(elementId2))
original_int_points = OS.ApplyInvTransform(ipTet)
target_int_points = TS.ApplyInvTransform(ipTet)
if self._debug_IntegrationMesh is not None: # pragma: no cover
AddElementToViz(self._debug_IntegrationMesh,original_int_points,ED.Point_1,"Left")
AddElementToViz(self._debug_IntegrationMesh,target_int_points,ED.Point_1,"Right")
# Append points and weights
weights_IPoints.extend(iwTet)
mesh_IElement.extend(np.full(len(tet_w),elementId1))
mesh_IIElement.extend(np.full(len(tet_w),elementId2))
for ip_nb in range(len(tet_w)):
meshI_IPoints.append(original_int_points[ip_nb,:])
meshII_IPoints.append(target_int_points[ip_nb,:])
continue
PrintProgressBar(len(selection1.indices),len(selection1.indices))
if self._debug_IntegrationMesh is not None :
self._debug_IntegrationMesh.PrepareForOutput()
self._debug_IntegrationMesh.GenerateManufacturedOriginalIDs()
weights_IPoints = np.array(weights_IPoints)
meshI_IPoints = np.array(meshI_IPoints)
meshII_IPoints= np.array(meshII_IPoints)
mesh_IElement = np.array(mesh_IElement)
mesh_IIElement = np.array(mesh_IIElement)
totalNumberOfIP = weights_IPoints.shape[0]
if len(meshI_IPoints) == 0:
print("Warning! -> Zero elements in contact")
return
# need to code the transfer of the field to the integration points meshI
from Muscat.MeshTools.MeshFieldOperations import GetFieldTransferOp
meshIOps = {}
# we must apply the offsets to the operators to build correctly the coupling terms.
for f, offset in zip(fields,offsetsI):
op, status, entities = GetFieldTransferOp(f,meshI_IPoints,method="Interp/Clamp",elementFilter=ElementFilter(eTag=self.on) )
# apply the offset to the op
op.col += offset
op.resize( (totalNumberOfIP,totalNumberOfDofsI) )
meshIOps[f.name] = (op, status)
# need to code the transfer of the field to the integration points meshII
meshIIOps = {}
for f,offset in zip(fieldsII,offsetsII):
op, status, entities = GetFieldTransferOp(f,meshII_IPoints,method="Interp/Clamp",elementFilter=ElementFilter(eTag=self.onII) )
op.col += offset
op.resize( (totalNumberOfIP,totalNumberOfDofsII) )
meshIIOps[f.name] = (op, status)
diagWMatrix = diags(weights_IPoints)
for arg in self.args:
fieldI = fieldDicI.get(arg)
fieldII = fieldDicII.get(arg)
opI = meshIOps[arg][0]
opII = meshIIOps[arg][0]
A = diagWMatrix.dot(opI).T.dot(opI)
B = diagWMatrix.dot(opI).T.dot(opII)
## different types of behaviors
if meshII is meshI:
## if 2 sides are in the same mesh (monolithic problem)
## the we add the constraint to the system
data = (A-B)
nums= np.where(np.sum(np.abs(data),axis=1)!= 0)[0]
CH.AddEquationsFromOperatorAAndb((data[nums,:]).tocoo())
else:
pass
## in this case we are in "advance mode"
## the user is responsible of using the operators
self.weights_IPoints = weights_IPoints
self.meshIOps = meshIOps
self.meshIIOps = meshIIOps
self.meshI_IPoints = meshI_IPoints
self.meshII_IPoints = meshII_IPoints
return CH
def __computeNormalSurface(self,nodes,conn):
barycenter = np.sum(nodes[conn ,:],axis=0)/len(conn)
edgeVector = nodes[conn[0],:] - nodes[conn[1],:]
planeVector = barycenter - nodes[conn[1],:]
normal = np.cross(edgeVector, planeVector)
normal /= np.linalg.norm(normal)
return normal, barycenter, planeVector/np.linalg.norm(planeVector)
def __computeNormalEdge(self,nodes,conn):
barycenter = np.sum(nodes[conn ,:],axis=0)/len(conn)
planeVector = barycenter - nodes[conn[0],:]
normal = np.array([planeVector[1], -planeVector[0],0])
normal /= np.linalg.norm(normal)
return normal, barycenter, planeVector/np.linalg.norm(planeVector)
def __ComputeNormal(self,subMesh,name,data,elId):
if ED.dimensionality[name] == 1:
return self.__computeNormalEdge(subMesh.nodes,data.connectivity[elId,:])
elif ED.dimensionality[name] == 2:
return self.__computeNormalSurface(subMesh.nodes,data.connectivity[elId,:])
else: # pragma: no cover
raise Exception(" Error ")
[docs]
def AreCCW(p1,p2,p3):
return bool((np.sign(np.cross(p2-p1,p3-p2)[2])+1)/2)
[docs]
def Append(l,point,tol):
p = np.array(point)
if len(l):
if np.linalg.norm(l[-1] - p ) > tol:
l.append(p)
else:
l.append(point)
[docs]
def CheckTriWinding(points, tri, allowReversed):
trisq = np.ones((3,3))
trisq[:,0:2] = np.array(points[tri,0:2])
detTri = np.linalg.det(trisq)
if detTri < 0.0:
if allowReversed:
#print(tri)
a = tri[2]
tri[2] = tri[1]
tri[1] = a
return tri
a = trisq[2,:].copy()
trisq[2,:] = trisq[1,:]
trisq[1,:] = a
else: # pragma: no cover
raise ValueError("triangle has wrong winding direction")
return tri
[docs]
def Intersection(points1, _poly1,points2,_poly2, tol):
## cut poly1 by poly2
poly1 = list(_poly1)
poly1 = CheckTriWinding(points1,poly1,True)
poly1 = [ points1[x,:] for x in poly1 ]
poly2 = list(_poly2)
poly2 = CheckTriWinding(points2,poly2,True)
poly2.append(_poly2[0])
res = []
for cuts in range(len(poly2)-1):
cp0 = points2[poly2[cuts]]
cp1 = points2[poly2[cuts+1]]
res = []
Append(poly1,poly1[0],tol)
for s in range(len(poly1)-1):
sp0 = poly1[s]
sp1 = poly1[s+1]
if AreCCW(cp0,cp1,sp0):
# point inside keep the point
Append(res,sp0,tol)
if AreCCW(cp0,cp1,sp0) != AreCCW(cp0,cp1,sp1):
# segment must be cut by cutter
# keep the intersection
x1 = cp0[0]
y1 = cp0[1]
x2 = cp1[0]
y2 = cp1[1]
x3 = sp0[0]
y3 = sp0[1]
x4 = sp1[0]
y4 = sp1[1]
if (( x1-x2 )*(y3-y4)-(y1-y2)*(x3-x4) ) == 0:
px = x4
py = y4
else:
px = ((x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4)) /(( x1-x2 )*(y3-y4)-(y1-y2)*(x3-x4) )
py = ((x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4)) /(( x1-x2 )*(y3-y4)-(y1-y2)*(x3-x4) )
px = np.clip(px, min(x3,x4),max(x3,x4))
py = np.clip(py, min(y3,y4),max(y3,y4))
Append(res,[px,py,0],tol)
if AreCCW(cp0,cp1,sp1):
# the last point is treated by the next integration
Append(res,sp1,tol)
poly1 = res
if len(poly1) < 3:
break
if len(res):
Append(res,res[0],tol)
return np.array(res)
[docs]
def IntersectionOf2ConvexHulls(pointsI,pointsII):
from scipy.spatial import ConvexHull, HalfspaceIntersection
from scipy.optimize import linprog
# computation of the convex Hulls
CHI = ConvexHull(points=pointsI)
CHII = ConvexHull(points=pointsII)
# computation of a interior_point
# using method from
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.HalfspaceIntersection.html
halfSpaces = np.vstack((CHI.equations, CHII.equations))
norm_vector = np.reshape(np.linalg.norm(halfSpaces[:, :-1], axis=1), (halfSpaces.shape[0], 1))
c = np.zeros((halfSpaces.shape[1],))
c[-1] = -1
A = np.hstack((halfSpaces[:, :-1], norm_vector))
b = - halfSpaces[:, -1:]
res = linprog(c, A_ub=A, b_ub=b, bounds=(None, None))
# no intersection return false
if res.success == False:
return False, None, None
x = res.x[:-1]
# computation of the intersection
try:
HI = HalfspaceIntersection(halfSpaces,interior_point = x, incremental=False)
except:
return False, None, None
# we use all the intersection point to build the intersected convex hull
CHIntersection = ConvexHull(points=HI.intersections)
nbPoints = HI.intersections.shape[0]
if nbPoints == 4:
""" we know the only possibility in 3D is a tetra """
points = HI.intersections
tets = np.array([np.hstack((CHIntersection.simplices[0],[x for x in range(4) if x not in CHIntersection.simplices[0]]))])
else:
nbPoints += 1
points = np.vstack((HI.intersections,np.mean(HI.intersections,axis=0)))
nbPoints = points.shape[0]
tets = np.empty((len(CHIntersection.simplices),4),dtype=int)
tets[:,0:3] = CHIntersection.simplices
tets[:,3] = nbPoints -1
## verification of the order of the tets
for i in range(tets.shape[0]):
p0 = points[tets[i,0],:]
p1 = points[tets[i,1],:]
p2 = points[tets[i,2],:]
p3 = points[tets[i,3],:]
if np.dot(p3-(p0+p1+p2)/3, np.cross(p1-p0,p2-p1)) < 0 :
# swap if the forth point is on the wrong side
tets[i,1], tets[i,0] = tets[i,0], tets[i,1]
return True, points, tets
[docs]
def CheckIntegrity3D(GUI:bool=False):
from Muscat.MeshTools.MeshCreationTools import CreateCube
nx,ny, nz = 1, 1, 1
meshI = CreateCube(dimensions=[nx+1,ny+1,nz+1],origin=[0,0,0], spacing=[10./nx,10./ny,1/nz], ofTetras=True)
meshII = CreateCube(dimensions=[nx+1,ny+1,nz+1],origin=[5,0,0], spacing=[1./nx,10./ny,10/nz], ofTetras=True)
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1
from Muscat.Actions.OpenInParaView import OpenInParaView
numberingI = ComputeDofNumbering(meshI,LagrangeSpaceGeo,fromConnectivity=True)
numberingII = ComputeDofNumbering(meshII,LagrangeSpaceGeo,fromConnectivity=True)
tag1 = "3D"
tag2 = "3D"
space = LagrangeSpaceP1
unknownNames = ["T"]
unknownNamesII = ["T"]
fieldsI = []
for name in unknownNames:
fieldsI.append(FEField(name,meshI,space,numberingI))
fieldsII = []
for name in unknownNamesII:
fieldsII.append(FEField(name,meshII,space,numberingII))
obj = KRMortar()
obj.AddArg("T")
obj.From()
obj.To(offset=[0,0,0] )
obj.SideI(tag1)
obj.SideII([tag2])
obj._debug_IntegrationMesh = Mesh()
print(obj)
print(obj.GenerateEquations(meshI,fieldsI,meshII=meshII,fieldsII=fieldsII))
if GUI: # pragma: no cover
from Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(mesh=meshI)
OpenInParaView(mesh=meshII)
OpenInParaView(mesh=obj._debug_IntegrationMesh)
return "ok"
[docs]
def CheckIntegrityIntersectionConvexHull3D(GUI:bool=False):
from scipy.spatial import ConvexHull
pointsI = np.array([[0.0,0.0,0.0],
[ 1.0,0.0,0.0],
[ 0.0,1.0,0.0],
[ 0.0,0.0,1.0]])
pointsII = np.array([[0.5,0.0,0.1],
[1.5,0,0],
[1.5,1,0],
[1.5,0,1]])
status, points, tets = IntersectionOf2ConvexHulls(pointsI,pointsII)
if status == False:# pragma: no cover
raise Exception("Error in the detection of the intersection")
mesh = Mesh()
mesh.nodes = np.vstack((pointsI,pointsII,points))
elements = mesh.GetElementsOfType(ED.Triangle_3)
CHI = ConvexHull(points=pointsI)
CHII = ConvexHull(points=pointsII)
CHintersection = ConvexHull(points=points)
cpt = 0
for simplex in CHI.simplices:
elements.AddNewElement(simplex,cpt)
cpt += 1
for simplex in CHII.simplices:
elements.AddNewElement(simplex+pointsI.shape[0],cpt)
cpt += 1
if tets.shape[0] != 1:# pragma: no cover
raise Exception("ERROR! Wrong number of tetras in the intersection ")
elements = mesh.GetElementsOfType(ED.Tetrahedron_4)
elements.connectivity = tets+pointsI.shape[0]+pointsII.shape[0]
elements.cpt = tets.shape[0]
mesh.GenerateManufacturedOriginalIDs()
mesh.PrepareForOutput()
if GUI:# pragma: no cover
from Muscat.Bridges.vtkBridge import PlotMesh
PlotMesh(mesh)
from Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(mesh=mesh)
return "ok"
[docs]
def CheckIntegrityIntersection(GUI:bool=False):
P = np.array([[0.,0,0],
[5.,0,0],
[0.,5,0],
[0,6,0]])
data = Intersection(P, [0,1,3],P,[0,1,2],tol=0.0001)
center = np.sum(data,axis=0)/(data.shape[0]-1)
result = (center,data)
P = np.array([[0.,0,0],
[4.,0,0],
[0.,3,0],
[2.,-1,0],
[3.,2,0],
[-1.,3,0]])
data = Intersection(P, [0,1,2],P,[3,4,5],tol=0.0000001)
center = np.sum(data[0:-1,:],axis=0)/(data.shape[0]-1)
result = (center,data)
return "ok"
[docs]
def CheckIntegrity1DInterface2Meshes(GUI:bool=False):
"""CheckIntegrity for 2 1D meshes
Mesh 1::
y ^
|
1|x 2 x
||
|| 3 x
||
0 -x-------4-x----> x
Mesh 2::
y ^
|
1 | 2 x
| |
| 3 x
| |
0 --------4-x----> x
"""
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1
ps = np.array([[0.,0.,0],
[0.,1.,0],
[1.,1,0],
[1.,0.5,0],
[1.,0.,0],
])
es = np.array([[0,1]])
meshI = CreateMeshOf(ps,es, ED.Bar_2)
ps = np.array([ [1.,1,0],
[1.,0.5,0],
[1.,0.,0],
])
es = np.array([[0,1],[1,2]])
meshII = CreateMeshOf(ps,es, ED.Bar_2)
tag1 = "Left"
tag2 = "Right"
meshI.GetElementsOfType(ED.Bar_2).tags.CreateTag(tag1).SetIds(0)
meshII.GetElementsOfType(ED.Bar_2).tags.CreateTag(tag2).SetIds([0,1])
print(meshI)
print(meshII)
numberingI = ComputeDofNumbering(meshI,LagrangeSpaceGeo)
numberingII = ComputeDofNumbering(meshII,LagrangeSpaceGeo)
print("nDofs I",numberingI.size)
print("nDofs II",numberingII.size)
print("------------------------")
space = LagrangeSpaceP1
unknownNames = ["T"]
fieldsI = []
fieldsII = []
for name in unknownNames:
fieldsI.append(FEField(name,meshI,space,numberingI))
fieldsII.append(FEField(name,meshII,space,numberingII))
obj = KRMortar()
obj.AddArg("T")
#obj.From()
obj.To(offset=[1,0,0] )
obj.SideI(tag1)
obj.SideII(tag2)
print(obj)
if GUI :
obj._debug_IntegrationMesh = Mesh()
CH = obj.GenerateEquations(meshI,fieldsI,CH=None,meshII=meshII,fieldsII=fieldsII)
if GUI : # pragma: no cover
from Muscat.Bridges.vtkBridge import PlotMesh
PlotMesh(meshI)
PlotMesh(meshII)
PlotMesh(obj._debug_IntegrationMesh)
return "ok"
[docs]
def CheckIntegrity1DInterface(GUI:bool=False):
"""
CheckIntegrity for::
y ^
|
1 x 2 x
| |
| 3 x
| |
0 -x-------4-x----> x
"""
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1
ps = np.array([[0.,0.,0],
[0.,1.,0],
[1.,1,0],
[1.,0.5,0],
[1.,0.,0],
])
es = np.array([[0,1],[2,3],[3,4]])
mesh = CreateMeshOf(ps,es, ED.Bar_2)
tag1 = "Left"
tag2 = "Right"
mesh.GetElementsOfType(ED.Bar_2).tags.CreateTag(tag1).SetIds(0)
mesh.GetElementsOfType(ED.Bar_2).tags.CreateTag(tag2).SetIds([1,2])
print(mesh)
numbering = ComputeDofNumbering(mesh,LagrangeSpaceGeo,fromConnectivity=True)
space = LagrangeSpaceP1
unknownNames = ["T"]
fields = []
for name in unknownNames:
fields.append(FEField(name,mesh,space,numbering))
obj = KRMortar()
obj.AddArg("T")
obj.To(offset=[1,0,0] )
obj.SideI(tag1)
obj.SideII(tag2)
obj._debug_IntegrationMesh = Mesh()
CH = obj.GenerateEquations(mesh,fields)
CH.SetNumberOfDofs(5)
print(CH.ToSparseFull().toarray())
if GUI : # pragma: no cover
from Muscat.Bridges.vtkBridge import PlotMesh
PlotMesh(mesh)
PlotMesh(obj._debug_IntegrationMesh)
return "ok"
[docs]
def CheckIntegrity2DScalar(GUI:bool=False):
from Muscat.TestData import GetTestDataPath
filename = GetTestDataPath()+"dent.msh"
from Muscat.IO.UniversalReader import ReadMesh
import Muscat.IO.GmshReader # to register the GmshReader
mesh = ReadMesh(filename)
#print(mesh.nodes.shape)
newNodes = np.zeros((mesh.nodes.shape[0],3),dtype=float)
newNodes[0:,:2] = mesh.nodes
mesh.nodes = newNodes
#print(mesh.nodes.shape)
#raise
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1
numbering = ComputeDofNumbering(mesh,LagrangeSpaceGeo,fromConnectivity=True)
tag1 = "Left"
tag2 = "Right"
space = LagrangeSpaceP1
unknownNames = ["T"]
fields = []
for name in unknownNames:
fields.append(FEField(name,mesh,space,numbering))
obj = KRMortar()
obj.AddArg("T")
obj.From(first=[-np.sin(np.pi/5),np.cos(np.pi/5),0 ], second=[-np.cos(np.pi/5),-np.sin(np.pi/5),0 ] )
angle = np.pi/5+0.01
obj.To(first=[np.sin(angle),np.cos(angle),0 ], second=[-np.cos(angle),np.sin(angle),0 ] )
obj.SideI(tag1)
obj.SideII(tag2)
obj._debug_IntegrationMesh = Mesh()
print(obj)
print(obj.GenerateEquations(mesh,fields))
if GUI: # pragma: no cover
from Muscat.Bridges.vtkBridge import PlotMesh
PlotMesh(obj._debug_IntegrationMesh)
return "ok"
[docs]
def CheckIntegrity3DVector(GUI:bool=False):
from Muscat.TestData import GetTestDataPath
filename = GetTestDataPath()+"dent3D.msh"
from Muscat.IO.UniversalReader import ReadMesh
import Muscat.IO.GmshReader
mesh = ReadMesh(filename)
#print(mesh)
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1
numbering = ComputeDofNumbering(mesh,LagrangeSpaceGeo,fromConnectivity=True)
tag1 = "Left"
tag2 = "Right"
space = LagrangeSpaceP1
unknownNames = ["u_0","u_1","u_2"]
fields = []
for name in unknownNames:
fields.append(FEField(name,mesh,space,numbering))
obj = KRMortar()
obj.AddArg("u_0")
obj.AddArg("u_1")
obj.AddArg("u_2")
obj.From(offset=[0,0,0], first=[-np.sin(np.pi/5),np.cos(np.pi/5),0 ], second=[-np.cos(np.pi/5),-np.sin(np.pi/5),0 ] )
angle = np.pi/5
obj.To (first=[np.sin(angle),np.cos(angle),0 ], second=[-np.cos(angle),np.sin(angle),0 ] )
obj.SideI(tag1)
obj.SideII(tag2)
obj._debug_IntegrationMesh = Mesh()
print(obj.GenerateEquations(mesh,fields))
if GUI: # pragma: no cover
from Muscat.Bridges.vtkBridge import PlotMesh
PlotMesh(obj._debug_IntegrationMesh)
return "ok"
[docs]
def CheckIntegrity(GUI:bool=False):
func = [CheckIntegrity1DInterface,
CheckIntegrity1DInterface2Meshes,
CheckIntegrity2DScalar,
CheckIntegrity3DVector,
CheckIntegrityIntersection,
CheckIntegrity3D,
CheckIntegrityIntersectionConvexHull3D
]
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities
return RunListOfCheckIntegrities(func, GUI=GUI)
if __name__ == "__main__": # pragma: no cover
print(CheckIntegrity(GUI=True))