# symfem-bilinearform-matrix-assembly.py
import meshzoo
import symfem
import numpy as np
from sympy import diff, symbols, integrate
def create_dof_map(points):
dof_map = {}
for i, point in enumerate(points):
dof_map[i] = point
return dof_map
def bilinear_form(grad_trial_x, grad_trial_y, grad_test_x, grad_test_y, κ):
return κ * (float(grad_trial_x)*float(grad_test_x) + float(grad_trial_y)*float(grad_test_y))
def global_K_assemble(points, cells, symfem_element):
x, y = symbols('x y')
global_stiffness_matrix = np.zeros((len(points), len(points)))
for cell in cells:
local_stiffness_matrix = np.zeros((len(cell), len(cell)))
map = points[cell]
local_basis_functions = symfem_element.map_to_cell(map)
##quadrature_points = [(1/3, 1/3)]
###quadrature_weights = [1/2]
quadrature_points = [(1/6, 1/6), (2/3, 1/6), (1/6, 2/3)]
quadrature_weights = [1/6, 1/6, 1/6]
area = abs((map[0][0]*(map[1][1]-map[2][1]) + map[1][0]*(map[2][1]-map[0][1]) + map[2][0]*(map[0][1]-map[1][1])))
for i, basis_i in enumerate(local_basis_functions):
for j, basis_j in enumerate(local_basis_functions):
integral = 0
for qp, w in zip(quadrature_points, quadrature_weights):
x,y=symbols('x y')
#grad_i = np.array([basis_i.diff(x).subs({x: qp[0], y: qp[1]}), basis_i.diff(y).subs({x: qp[0], y: qp[1]})])
grad_i = np.array([basis_i.diff(x).subs('x', qp[0]).subs('y', qp[1]), basis_i.diff(y).subs('x', qp[0]).subs('y', qp[1])])
#grad_j = np.array([basis_j.diff(x).subs({x: qp[0], y: qp[1]}), basis_j.diff(y).subs({x: qp[0], y: qp[1]})])
grad_j = np.array([basis_j.diff(x).subs('x', qp[0]).subs('y', qp[1]), basis_j.diff(y).subs('x', qp[0]).subs('y', qp[1])])
grad_product = np.dot(grad_i, grad_j)
integral += float(w) * float(grad_product)
local_stiffness_matrix[i, j] += integral * area
for i, global_i in enumerate(cell):
for j, global_j in enumerate(cell):
global_stiffness_matrix[global_i, global_j] += local_stiffness_matrix[i, j]
return global_stiffness_matrix
points, cells = meshzoo.rectangle_tri(
np.linspace(0.0, 4.0, 3),
np.linspace(0.0, 4.0, 3),
variant="zigzag"
)
dof_map = create_dof_map(points)
for dof, point in dof_map.items():
print(f"Degree of freedom {dof} corresponds to point {point}")
κ = 0.83
p1 = symfem.create_element("triangle", "Lagrange", 1)
Pn = p1.get_polynomial_basis()
print('Pn:', Pn)
gradient_list = []
for cell in cells:
map = points[cell]
cell_basis_functions = p1.map_to_cell(map)
cell_grads = [[diff(bf, x) for bf in cell_basis_functions] for x in symbols('x y')]
gradient_list.append(cell_grads)
κ = 0.83
bilinear_form_list = np.zeros((len(points), len(points)))
quadrature_points = [(1/6, 1/6), (2/3, 1/6), (1/6, 2/3)]
quadrature_weights = [1/6, 1/6, 1/6]
for cell_index, cell in enumerate(cells):
map = points[cell]
area = abs((map[0][0]*(map[1][1]-map[2][1]) + map[1][0]*(map[2][1]-map[0][1]) + map[2][0]*(map[0][1]-map[1][1])))
for i in range(len(cell)):
for j in range(len(cell)):
integral = 0
for qp, w in zip(quadrature_points, quadrature_weights):
grad_trial_x = gradient_list[cell_index][0][i].subs('x', qp[0]).subs('y', qp[1])
grad_trial_y = gradient_list[cell_index][1][i].subs('x', qp[0]).subs('y', qp[1])
grad_test_x = gradient_list[cell_index][0][j].subs('x', qp[0]).subs('y', qp[1])
grad_test_y = gradient_list[cell_index][1][j].subs('x', qp[0]).subs('y', qp[1])
integral += w * bilinear_form(grad_trial_x, grad_trial_y, grad_test_x, grad_test_y, κ)
bilinear_form_list[cell[i], cell[j]] += integral * area
print('bilinear_form_matrix:', bilinear_form_list)
K = global_K_assemble(points, cells, p1)
print('K:', K)
#def load_function(x, y):
#return 10
def jacobian_det(map):
x1, y1 = map[0]
x2, y2 = map[1]
x3, y3 = map[2]
return abs(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))
global_linear_form_vector = np.zeros(len(points))
gaussian_weights = [1/6, 1/6, 1/6] # weights for 3-point Gaussian quadrature on a triangle
for cell in cells:
local_linear_form_vector = np.zeros(len(cell))
locations = points[cell]
detJ = jacobian_det(locations)
mapped_points = p1.map_to_cell(locations)
for i in range(3):
mpf=mapped_points[i].subs('x', locations[i][0]).subs('y', locations[i][1])
integral = 10*float(mpf) * float(gaussian_weights[i]) * float(detJ)
local_linear_form_vector[i] = integral
global_linear_form_vector[cell] += local_linear_form_vector
print('global_linear_form_vector:', global_linear_form_vector)
rhs = K * global_linear_form_vector
solution = np.linalg.solve(bilinear_form_list, rhs)
print('solution:')
print()
print('length:', len(solution))
print('length_row:', len(solution[0]))
np.set_printoptions(precision=4)
print(solution)
print()
def map_solution_to_nodes(solution, dof_map):
"""
Maps the solution vector to the nodes of the rectangle mesh.
Parameters:
solution (numpy.ndarray): The solution vector of size 9x1.
dof_map (dict): The degree of freedom map that maps each node in the mesh to a unique integer index.
Returns:
dict: A dictionary that maps each node in the mesh to its corresponding solution value.
"""
node_solutions = { }
for dof, point in dof_map.items():
node_solutions[tuple(point)] = solution[dof]
return node_solutions
node_solutions = map_solution_to_nodes(solution, dof_map)
for node, solution_value in node_solutions.items():
#solution_value_formatted = "{:.4f}".format(solution_value)
solution_value_formatted = [ "{:.4f}".format(val) for val in solution_value ]
print(f"Node {node} corresponds to solution value {solution_value_formatted}")
So far what I feel unsure about is how to correctly set a starting temperature and also how to make sense out of the solution. The mesh contains 9 nodes and 8 cells and so far as I can tell the solution is a 9x9 list (or matrix if so to call it.). So far my intent is to apply a force of 10 joules at each node. Anyone can assist me to complete such an example for steady state heat transfer?