FW: [EXT] Re: rigid contact with arbitrary surface

Actually I have a python code for distance to surface

Here I have about 15 points in my brain placed on mri and I want to know the 3 ones closer to the cortex

#############################################################

#!/venv/bin/python3

# -*- coding: UTF8 -*-

# COPYRIGHT

# MD Anderson Cancer Center

# Imaging Physics department

# Last modified 03/03/22

# Anne-Cecile Lesage. Research scientist

import sys

import numpy as np

import math as ma

from numpy import linalg as LA

# The algorithm is based on

# math.stackexchange.com/questions/588871/minimum-distance-between-point-and-face

def mindistptri(p,p1,p2,p3,n,t,p0,dv,a,b):

x1 = p1[0]

y1 = p1[1]

z1 = p1[2]

x2 = p2[0]

y2 = p2[1]

z2 = p2[2]

x3 = p3[0]

y3 = p3[1]

z3 = p3[2]

#nx = (y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)

#ny = (z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)

#nz = (x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)

#norm = np.sqrt(nx*nx+ny*ny+nz*nz)

#n[0] = nx/norm

#n[1] = ny/norm

#n[2] = nz/norm

t[:] = n[:]*p1[:]-n[:]*p[:]

p0[:] = p[:]+n[:]*t[:]

dmin = 1000000

dv[:] = p[:]-p0[:]

dp = LA.norm(dv[:])

a[0:3,0]=p1[0:3]

a[0:3,1]=p2[0:3]

a[0:3,2]=p3[0:3]

b[:] = p0[:]

x = np.linalg.solve(a,b)

dmin=1000000

d=dmin

compt = 0

if x[0]>0.0 and x[0]<1.0:

compt = compt + 1

# print(x)

pointin = 0

if compt==3:

pointin = 1

d = dp

print("Point projection inside triangle\n")

elif compt<3:

dv[:] = p[:]-p1[:]

dp = LA.norm(dv[:])

if dp<dmin:

dmin=dp

dv[:] = p[:]-p2[:]

dp = LA.norm(dv[:])

if dp<dmin:

dmin=dp

dv[:] = p[:]-p3[:]

dp = LA.norm(dv[:])

if dp<dmin:

dmin=dp

d = dmin

#print(compt)

#print(pointin)

return(d)

#################################################################

print('Read 3 pois parameter file')

fileSim = "3pois.dat"

fileObj = open(fileSim)

params = {}

for line in fileObj:

line = line.strip()

#print(line)

keyvalue = line.split("=")

#print(keyvalue)

if len(keyvalue) == 2:

params[keyvalue[0].strip()] = keyvalue[1].strip()

fileObj.close()

# print dictionnary

#print(params)

# Convert the read parameters into desired types

landmarkdata = params["LANDMARKDATA"]

brainsestlfile = params["BRAINSESTLFILE"]

distfile = params["DISTFILE"]

print(landmarkdata,' ',brainsestlfile,' ',distfile,'\n')

####################################################################################

print("Part 1 read landmarks\n")

filename = "%s.txt" % landmarkdata

file = open(filename,"r")

case25 = np.loadtxt(file,dtype=float)

file.close()

nl = np.size(case25,0)

#print(case25)

preopt = np.zeros((nl,3))

preopt[:,0:2] = case25[:,0:2]*-0.001

preopt[:,2] = case25[:,2]*0.001

#####################################################################################################################################

# read stl file:

print('Read stl file brain exterior surface self written subroutine\n')

filename = "%s.stl" % brainsestlfile

#print(filename)

fileObj = open(filename)

ntri = 0

for line in fileObj:

line = line.strip()

#print(line)

keyvalue = line.split(" ")

#print(keyvalue)

if len(keyvalue) == 5:

#print(line)

#print(keyvalue)

ntri = ntri + 1

fileObj.close()

ptstl = np.zeros((ntri*3,3))

normstl = np.zeros((ntri,3))

normvstlr = np.zeros((ntri,3))

fileObj = open(filename)

ntri = 0

npstl = 0

for line in fileObj:

line = line.strip()

#print(line)

keyvalue = line.split(" ")

#print(keyvalue)

if len(keyvalue) == 5:

normstl[ntri,0] = float(keyvalue[2])

normstl[ntri,1] = float(keyvalue[3])

normstl[ntri,2] = float(keyvalue[4])

#print(normstl[ntri,:])

ntri = ntri + 1

#print(keyvalue)

if len(keyvalue) == 4:

ptstl[npstl,0] = float(keyvalue[1])

ptstl[npstl,1] = float(keyvalue[2])

ptstl[npstl,2] = float(keyvalue[3])

npstl = npstl + 1

fileObj.close()

ptstl[:,:] = ptstl[:,:]*0.001

print(" ntri ", ntri," npstl ", npstl,"\n")

#norm = np.zeros(ntri)

#for j in range(10):

# norm = LA.norm(normstl[j,:])

# print(" tri ",j+1," norm ", norm,"\n")

# The algorithm is based on

# math.stackexchange.com/questions/588871/minimum-distance-between-point-and-face

# toy example

p = np.zeros(3)

p1 = np.zeros(3)

p2 = np.zeros(3)

p3 = np.zeros(3)

n = np.zeros(3)

t = np.zeros(3)

p0 = np.zeros(3)

dv = np.zeros(3)

a = np.zeros((3,3))

b = np.zeros(3)

x = np.zeros(3)

dist = np.zeros(nl)

p1[0] = 1.0

p1[1] = 1.0

p1[2] = 1.0

p2[0] = 2.0

p2[1] = 1.0

p2[2] = 1.0

p3[0] = 1.0

p3[1] = 2.0

p3[2] = 1.0

p[0] = 20.0

p[1] = 11.0

p[2] = 0.5

#p[0] = 1.1

#p[1] = 1.1

#p[2] = 0.5

dmin = 100000

filedist = '%s.dat' % distfile

dmin3l1 =100000

l1 = 0

l2 = 0

l3 = 0

with open(filedist, "w") as file1:

for i in range(nl):

p[0:3] = preopt[i,0:3]

npstl = 0

dmin = 1000000

for j in range(ntri):

p1[0:3]= ptstl[npstl,0:3]

npstl = npstl + 1

p2[0:3]= ptstl[npstl,0:3]

npstl = npstl + 1

p3[0:3]= ptstl[npstl,0:3]

npstl = npstl + 1

n[0] = normstl[j,0]

n[1] = normstl[j,1]

n[2] = normstl[j,2]

d=mindistptri(p,p1,p2,p3,n,t,p0,dv,a,b)

if d<dmin:

dmin = d

dist[i] = dmin

if dmin < dmin3l1:

dmin3l1 = dmin

l1 = i+1

file1.write("%d %.6f\n" % (i+1, dmin))

print("landmark ",i+1," d ",dmin)

dmin3l2 =100000

for i in range(nl):

if dist[i] < dmin3l2 and dist[i]> dmin3l1:

dmin3l2 = dist[i]

l2 = i+1

dmin3l3 =100000

for i in range(nl):

if dist[i] < dmin3l3 and dist[i]> dmin3l2:

dmin3l3 = dist[i]

l3 = i+1

print(" l1 ",l1," l2 ",l2," l3 ",l3,"\n")

From: Lesage,Anne Cecile J
Sent: Monday, March 28, 2022 6:37 PM
To: Konstantinos Poulios <logari81@googlemail.com>
Subject: RE: [EXT] Re: rigid contact with arbitrary surface

Dear Konstantinos

I shall be able to generate the obs array exteriorly

What do you mean by signed distance?

I often simulate organs that can move inside a 2D surface cavity

So the distance shall have always the same sign isn’t it?

I have worked a lot with level sets.

So the algorithm distance from a point to a surface (surface being defined by a triangle collection) I know it almost by heart

I could implement it in getfem

I am aiming to study the source code of getfem these next weeks

My principal investigator is requiring that I develop the competence to change getfem

So that there is no limit to what I can do

Just that I need to refresh my C++ with tutorials

Rigid contact is a way to start easy with mastering contact methods in getfem

However my ultimate goal is to be able to model contact between two 3D objects but I did not have any success with it

My python script lines are those 3 options

Option 1 generates a core file

Option 2 is not the correct syntax

Option 3 takes forever to converge and complains about tetrahedra orientation

optslcont=3

print(" optslcont ", optslcont)

if optslcont==1:

md.add_penalized_contact_between_nonmatching_meshes_brick(mimu3h, "uh", "ub", "penal_r", HEAD_BOUND, BRAIN_BOUND)

elif optslcont==2:

md.add_integral_contact_between_nonmatching_meshes_brick(mimu3h, "uh", "ub", "penal_r", HEAD_BOUND, BRAIN_BOUND)

elif optslcont==3:

release_dist = 1.

aug_factor = 0.1;

md.add_initialized_data('r', aug_factor * Muh * (3*Lambdah + 2*Muh) / (Lambdah + Muh) )

md.add_initialized_data('f_coeff', 0.)

md.add_initialized_data('alpha', 1.)

ibc = md.add_integral_large_sliding_contact_brick_raytracing('r', release_dist, 'f_coeff', 'alpha', 0)

md.add_filtered_fem_variable('multcontactb', mfub, BRAIN_BOUND)

md.add_slave_contact_boundary_to_large_sliding_contact_brick(ibc, mimu3b, BRAIN_BOUND, 'ub', 'multcontactb', '')

md.add_master_contact_boundary_to_large_sliding_contact_brick(ibc, mimu3h, HEAD_BOUND, 'uh', '')

Regards

Thank you

Anne-Cecile

From: Konstantinos Poulios <logari81@googlemail.com>
Sent: Monday, March 28, 2022 5:27 PM
To: Lesage,Anne Cecile J <AJLesage@mdanderson.org>
Cc: getfem-users@nongnu.org
Subject: [EXT] Re: rigid contact with arbitrary surface

WARNING: This email originated from outside of MD Anderson. Please validate the sender's email address before clicking on links or attachments as they may not be safe.

Dear Anne-Cecile

A good starting point for modelling contact in GetFEM is this example:

https://git.savannah.nongnu.org/cgit/getfem.git/tree/interface/tests/python/demo_static_contact.py

The special thing in your case is that you would like to have contact between a 3D solid and a 2D surface, which AFAIR is a case that is not directly available, but everything is possible.

For example if you are interested in small-deformation (i.e. linearized kinematics) then this is how contact is defined, in the demo I gave you above

OBS = mfd.eval(obstacle)
md.add_initialized_fem_data('obstacle', mfd, OBS)
md.add_integral_contact_with_rigid_obstacle_brick(mim_friction, 'u', 'lambda_n',
'obstacle', 'r', GAMMAC, version-4);

In your case, "u" will be a displacement variable describing the deformation of the solid, "lambda_n" will be a scalar multiplier defined on the surface of the (deformable) 3D solid, mfd will be a scalar mesh_fem, defined on the same (3D solid) mesh as "u".

The only special thing for your case will be how you calculate the array OBS. This array must contain the (signed) distance of every node of (the 3D solid) mesh_fem mfd from your rigid 2D obstacle surface. AFAIR GetFEM does not have a convenience function for calculating the distance between points in a 3D solid mesh and a 2D surface mesh. The raytracing interpolation available in GetFEM allows only to calculate distances from the surface of a 3D solid to the surface of another 3D solid. So unless you find another way to calculate your OBS array, we will need to do some modification in GetFEM to make this easier.

Best regards

Kostas

On Mon, Mar 28, 2022 at 10:48 PM Lesage,Anne Cecile J <AJLesage@mdanderson.org> wrote:

Dear all

Is it possible in getfem to impose a contact with a rigid arbitrary surface defined by a 2D triangle mesh for example?

Regards

Anne-Cecile

The information contained in this e-mail message may be privileged, confidential, and/or protected from disclosure. This e-mail message may contain protected health information (PHI); dissemination of PHI should comply with applicable federal and state laws. If you are not the intended recipient, or an authorized representative of the intended recipient, any further review, disclosure, use, dissemination, distribution, or copying of this message or any attachment (or the information contained therein) is strictly prohibited. If you think that you have received this e-mail message in error, please notify the sender by return e-mail and delete all references to it and its contents from your systems.

From:	Lesage,Anne Cecile J
Subject:	FW: [EXT] Re: rigid contact with arbitrary surface
Date:	Tue, 29 Mar 2022 00:34:02 +0000