|
From: | Lesage,Anne Cecile J |
Subject: | FW: [EXT] Re: rigid contact with arbitrary surface |
Date: | Tue, 29 Mar 2022 00:34:02 +0000 |
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 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>
Dear Anne-Cecile A good starting point for modelling contact in GetFEM is this example: 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) 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:
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.
|
[Prev in Thread] | Current Thread | [Next in Thread] |