Skip to content

Commit

Permalink
fix: modifying for 3d/generic case.
Browse files Browse the repository at this point in the history 
2d supports are going to need to be updated
  • Loading branch information
Lachlan Grose committed Feb 14, 2022
1 parent a5d27f1 commit c6a8ea3
Showing 1 changed file with 28 additions and 85 deletions.
113 changes: 28 additions & 85 deletions LoopStructural/interpolators/_p1interpolator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,6 @@ def __init__(self, mesh):
self.shape = "rectangular"
DiscreteInterpolator.__init__(self, mesh)
# whether to assemble a rectangular matrix or a square matrix
self.interpolator_type = "P1"
self.nx = len(self.support.nodes[self.region])
self.support = mesh

self.interpolation_weights = {
Expand All @@ -48,14 +46,13 @@ def add_gradient_ctr_pts(self, w=1.0):
def add_norm_ctr_pts(self, w=1.0):
points = self.get_norm_constraints()
if points.shape[0] > 0:
grad, elements = self.support.evaluate_shape_derivatives(points[:, :2])
inside = elements > -1
area = self.support.element_area(elements[inside])
wt = np.ones(area.shape[0])
wt *= w * area
grad, elements, inside = self.support.evaluate_shape_derivatives(points[:, :3])
size = self.support.element_size[inside]
wt = np.ones(size.shape[0])
wt *= w * size
# print(grad[inside,:,:].shape)
# print(self.support.elements[elements[inside]].shape)
elements = np.tile(self.support.elements[elements[inside]], (2, 1, 1))
elements = np.tile(self.support.elements[elements[inside]], (3, 1, 1))

elements = elements.swapaxes(0, 1)
# elements = elements.swapaxes(0,2)
Expand All @@ -73,50 +70,43 @@ def add_norm_ctr_pts(self, w=1.0):
def add_ctr_pts(self, w=1.0):
points = self.get_value_constraints()
if points.shape[0] > 1:
N, tri = self.support.evaluate_shape(points[:, :2])
mask = tri > 0
area = self.support.element_area(tri[mask])
wt = np.ones(area.shape[0])
wt *= w * area
N, elements, inside = self.support.evaluate_shape(points[:, :3])
size = self.support.element_size[elements[inside]]
wt = np.ones(size.shape[0])
wt *= w * size
self.add_constraints_to_least_squares(
N[mask, :] * wt[:, None],
points[mask, 3] * wt,
self.support.elements[tri[mask], :],
N[inside, :] * wt[:, None],
points[inside, 3] * wt,
self.support.elements[elements[inside], :],
name="value",
)

def minimize_edge_jumps(
self, stren, w=0.1, maxmdDist=None, vector_func=None
self, w=0.1, vector_func=None
): # NOTE: imposes \phi_T1(xi)-\phi_T2(xi) dot n =0
# iterate over all triangles
# flag inidicate which triangles have had all their relationships added
v1 = self.support.nodes[self.support.edges][:, 0, :]
v2 = self.support.nodes[self.support.edges][:, 1, :]
ncp = 2
# cp = np.zeros((v1.shape[0],ncp,2))
# cp[:,0] = 0.25*v1 + 0.75*v2
# cp[:,1] = 0.27*v1 + 0.25*v2
bc_t1 = self.support.barycentre(self.support.edge_relationships[:, 0])
bc_t2 = self.support.barycentre(self.support.edge_relationships[:, 1])

v = v1 - v2
e_len = np.linalg.norm(v, axis=1)
normal = np.array([v[:, 1], -v[:, 0]]).T
normal /= np.linalg.norm(normal, axis=1)[:, None]
v1 = self.support.nodes[self.support.shared_elements][:, 0, :]
v2 = self.support.nodes[self.support.shared_elements][:, 1, :]
bc_t1 = self.support.barycentre[self.support.shared_element_relationships[:, 0]]
bc_t2 = self.support.barycentre[self.support.shared_element_relationships[:, 1]]
norm = self.support.shared_element_norm
shared_element_size = self.support.shared_element_size

# evaluate normal if using vector func for cp2
if vector_func:
normal = vector_func((v1 + v2) / 2)
norm = vector_func((v1 + v2) / 2)
# evaluate the shape function for the edges for each neighbouring triangle
Dt, tri1 = self.support.evaluate_shape_derivatives(
bc_t1, elements=self.support.edge_relationships[:, 0]
bc_t1, elements=self.support.shared_element_relationships[:, 0]
)
Dn, tri2 = self.support.evaluate_shape_derivatives(
bc_t2, elements=self.support.edge_relationships[:, 1]
bc_t2, elements=self.support.shared_element_relationships[:, 1]
)

# constraint for each cp is triangle - neighbour create a Nx12 matrix
const_t = np.einsum("ij,ikj->ik", normal, Dt)
const_n = -np.einsum("ij,ikj->ik", normal, Dn)
const_t = np.einsum("ij,ijk->ik", norm, Dt)
const_n = -np.einsum("ij,ijk->ik", norm, Dn)
# const_t_cp2 = np.einsum('ij,ikj->ik',normal,cp2_Dt)
# const_n_cp2 = -np.einsum('ij,ikj->ik',normal,cp2_Dn)

Expand All @@ -127,58 +117,11 @@ def minimize_edge_jumps(
# tri_cp2 = np.hstack([self.support.elements[cp2_tri1],self.support.elements[tri2]])
# add cp1 and cp2 to the least squares system
self.add_constraints_to_least_squares(
const * e_len[:, None] * w,
const * shared_element_size[:, None] * w,
np.zeros(const.shape[0]),
tri_cp1,
name="edge jump cp1",
name="edge jump",
)
# p2.add_constraints_to_least_squares(const_cp2*e_len[:,None]*w,np.zeros(const_cp1.shape[0]),tri_cp2, name='edge jump cp2')

def minimize_gradient_norm(
self, stren, w=0.1, maxmdDist=None
): # NOTE: imposes \phi_T1(xi)-\phi_T2(xi) dot n =0
# iterate over all triangles
# flag inidicate which triangles have had all their relationships added
v1 = self.support.nodes[self.support.edges][:, 0, :]
v2 = self.support.nodes[self.support.edges][:, 1, :]
ncp = 2
# cp = np.zeros((v1.shape[0],ncp,2))
# cp[:,0] = 0.25*v1 + 0.75*v2
# cp[:,1] = 0.27*v1 + 0.25*v2
bc_t1 = self.support.barycentre(self.support.edge_relationships[:, 0])
bc_t2 = self.support.barycentre(self.support.edge_relationships[:, 1])

v = v1 - v2
e_len = np.linalg.norm(v, axis=1)
normal = np.array([v[:, 1], -v[:, 0]]).T
normal /= np.linalg.norm(normal, axis=1)[:, None]
# evaluate the shape function for the edges for each neighbouring triangle
Dt, tri1 = self.support.evaluate_shape_derivatives(
bc_t1, elements=self.support.edge_relationships[:, 0]
)
Dn, tri2 = self.support.evaluate_shape_derivatives(
bc_t2, elements=self.support.edge_relationships[:, 1]
)

# constraint for each cp is triangle - neighbour create a Nx12 matrix

# const_t = np.einsum('ij,ikj->ik',normal,Dt)
# const_n = -np.einsum('ij,ikj->ik',normal,Dn)
# const_t_cp2 = np.einsum('ij,ikj->ik',normal,cp2_Dt)
# const_n_cp2 = -np.einsum('ij,ikj->ik',normal,cp2_Dn)

const = np.hstack([Dt, -Dn]).T
const = const.swapaxes(1, 2)

# get vertex indexes
tri_cp1 = np.hstack([self.support.elements[tri1], self.support.elements[tri2]])
tri_cp1 = np.tile(tri_cp1, (2, 1, 1))
B = np.zeros((const.shape[0], const.shape[1]))
# np.tile(tri_cp1,(2,1,1))
print(tri_cp1.shape, const.shape, B.shape)
# tri_cp2 = np.hstack([self.support.elements[cp2_tri1],self.support.elements[tri2]])
# add cp1 and cp2 to the least squares system
self.add_constraints_to_least_squares(
const * e_len[None, :, None] * w, B, tri_cp1, name="constant norm"
)
# p2.add_constraints_to_least_squares(const_cp2*e_len[:,None]*w,np.zeros(const_cp1.shape[0]),tri_cp2, name='edge jump cp2')

0 comments on commit c6a8ea3

Please sign in to comment.