Things are a little disheveled there in my last reply. This is a much more adequate, Step 1… to developing the scenario I have in mind. So this will run and attempt so far to build the matrix:
$$ \begin{bmatrix} \text{NH}{dx} & 0 & 0 & \text{RH}{dx} & 0 & 0 \\ 0 & \text{NH}{dy} & 0 & 0 & \text{RH}{dy} & 0 \\ 0 & 0 & \text{NH}{dz} & 0 & 0 & \text{NH}{dz} \\ \text{NH}{dy} & \text{NH}{dx} & 0 & \text{RH}{dy} & \text{RH}{dx} & 0 \\ 0 & \text{NH}{dz} & \text{NH}{dy} & 0 & \text{RH}{dz} & \text{RH}{dy} \\ \text{NH}{dz} & 0 & \text{NH}{dx} & \text{RH}{dz} & 0 & \text{RH}{dx} \end{bmatrix} $$
Is there a more polite way of saying this in dolfinx speak by using the dx operater maybe? So far if this can accomplish the goal for the mixed element situation I would like to set up a mock-up to solve this. The boundary conditions and the rest of the boiler plate can be taken from the templates I linked to… So far this is the best that I could to…
import dolfinx
from dolfinx import mesh, fem, plot
from mpi4py import MPI
import basix
import numpy as np
from petsc4py import PETSc
import ufl
from dolfinx.fem.petsc import LinearProblem
from dolfinx.fem import dirichletbc, locate_dofs_topological
from ufl import Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, inner, dot
import pyvista as pv
import sympy as sp
# Material property matrix for wood
def D_matrix(mu, lamb):
D = sp.Matrix([
[2*mu+lamb, lamb, lamb, 0, 0, 0],
[lamb, 2*mu+lamb, lamb, 0, 0, 0],
[lamb, lamb, 2*mu+lamb, 0, 0, 0],
[0, 0, 0, mu, 0, 0],
[0, 0, 0, 0, mu, 0],
[0, 0, 0, 0, 0, mu]
])
return D
# Material properties
E = 2000000.0 # Young's modulus in PSI (for structural wood)
nu = 0.299999 # Poisson's ratio (for structural wood)
L = 72 # Beam length in inches
W = 4.0 # Beam width in inches
w = 100.0 # Distributed load in lbs/inch
# Define the mesh with tetrahedral cells
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
[20, 6, 6], cell_type=mesh.CellType.tetrahedron)
# Define elements for displacements (3 DOF) and rotations (3 DOF)
elem_disp = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(3,))
elem_rot = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(3,))
elem_mixed = basix.ufl.mixed_element([elem_disp, elem_rot])
V = fem.functionspace(domain, elem_mixed)
# Not in use... Strain tensor as a symmetric gradient of the displacement field
def epsilon(u):
return 0.5 * (ufl.grad(u) + ufl.grad(u).T)
# Not in use... Stress tensor derived from Hooke's law
def sigma(u):
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu)) # First Lamé parameter
mu = E / (2 * (1 + nu)) # Shear modulus (second Lamé parameter)
return lambda_ * ufl.tr(epsilon(u)) * ufl.Identity(3) + 2 * mu * epsilon(u)
# Define trial and test functions for displacements (u) and rotations (theta)
(u, theta) = TrialFunctions(V)
(v, tau) = TestFunctions(V)
# Apply boundary conditions and identify subspaces
V0 = V.sub(0) # Displacement space
V1 = V.sub(1) # Rotation space
Q0, _ = V0.collapse()
Q1, _ = V1.collapse()
fdim = domain.topology.dim - 1
# Define the trial and test functions for the displacement part (translations only)
du = ufl.TrialFunction(V0) # Trial function for translations
u_ = ufl.TestFunction(V0) # Test function for translations
# Define the bilinear form for translations (displacement part)
a_translations = ufl.inner(ufl.grad(du), ufl.grad(u_)) * ufl.dx
# Assemble the bilinear form for translations
A_translations = fem.petsc.assemble_matrix(fem.form(a_translations))
A_translations.assemble()
# Define the trial and test functions for rotations (theta)
dtheta = ufl.TrialFunction(V1) # Trial function for rotations
theta_ = ufl.TestFunction(V1) # Test function for rotations
# Define the bilinear form for rotations (using gradient of rotation)
a_rotations = ufl.inner(ufl.grad(dtheta), ufl.grad(theta_)) * ufl.dx
# Assemble the bilinear form for rotations
A_rotations = fem.petsc.assemble_matrix(fem.form(a_rotations))
A_rotations.assemble()
# Print matrix sizes
print(f"Matrix size for translations: {A_translations.size[0]} x {A_translations.size[1]}")
print(f"Matrix size for rotations: {A_rotations.size[0]} x {A_rotations.size[1]}")
# Assuming A_rotations is already defined as a PETSc matrix
#attributes = [attr for attr in dir(A_rotations) if not attr.startswith('_')]
#print(attributes)
# Define the total number of DOF (degrees of freedom)
total_dof = 6174
# Create empty matrix of the required size (6174x6174)
joined_components_K = np.zeros((total_dof, total_dof))
# You are looping through the matrix to populate the 3x3 blocks
# Iterate over blocks as you did previously, but now update the larger matrix
for i in range(0, A_translations.size[0] - 12, 12):
start_i = i
end_i = start_i + 12
trans_block = A_translations[start_i:end_i, start_i:end_i]
rot_block = A_rotations[start_i:end_i, start_i:end_i]
# Loop through the diagonal blocks (each 3x3) for both translation and rotation
for j in range(0, 12, 3): # j will be the starting index for each 3x3 block
# Extract the 3x3 blocks for translation and rotation
NH_block = trans_block[j:j+3, j:j+3] # 3x3 translation block
RH_block = rot_block[j:j+3, j:j+3] # 3x3 rotation block
# Extract the individual components from NH_block
NHdx = NH_block[0, 0] # Translation in x-direction (Node 1)
NHdy = NH_block[1, 1] # Translation in y-direction (Node 1)
NHdz = NH_block[2, 2] # Translation in z-direction (Node 1)
# Extract the individual components from RH_block
RHdx = RH_block[0, 0] # Rotation about x-axis (Node 1)
RHdy = RH_block[1, 1] # Rotation about y-axis (Node 1)
RHdz = RH_block[2, 2] # Rotation about z-axis (Node 1)
# Create the 6x6 matrix from the blocks (NH and RH)
joined_components_Mat = np.array([
[NHdx, 0, 0, RHdx, 0, 0],
[0, NHdy, 0, 0, RHdy, 0],
[0, 0, NHdz, 0, 0, NHdz],
[NHdy, NHdx, 0, RHdy, RHdx, 0],
[0, NHdz, NHdy, 0, RHdz, RHdy],
[NHdz, 0, NHdx, RHdz, 0, RHdx]
])
# Calculate the global indices (map the local 6x6 block to the larger matrix)
global_i = start_i + j
global_j = start_i + j
# Insert the 6x6 block into the larger matrix (block_matrix_6174x6174)
# I guess this should really be called B, since it is a network of B's
# that can be looped over in the fashion of K building in the previous
# posts, where the transpose of B is used to multiply by D and the det_J....
joined_components_K[global_i:global_i+6, global_j:global_j+6] = joined_components_Mat
# Now you have a block matrix of dimensions 6174x6174
# You can print or process the matrix as needed
print(joined_components_K.shape)
So essentially Ku=F with 6 dofs per node can now be set up manually either by converting to NumPy or some other acceptible way. I would guess so far that there is a way to have dolfinx apply boundary conditions to the bi-linear form to make that more simple. Mostly how to set this up with UFL and how to solve it without doing all of this so verbosely in the manual sense is a source of confusion for me still…
So the B matrix network now needs to be looped over in the fashion mentioned in its comment. So there will be different types of elements and also Guassian schemes… Is there some facet or function of dolfinx or UFL, that can help to ease the pain and over-verbosity of doing this"?
So the joined K is now a system of B matrices… So about the use of the Jacobian there… When dolfinx builds the bi-linear form does it already apply the determinant of the Jacobian or does that have to be done manually? So far I take things as that the bi-linear form that dolfinx builds is the gradients of the shape functions evaluated at positions which canverts to a scalar in local space? How does that differ from the gradients of the shape functions evaluated at positions in which converts into a scalar multiplied by the determinant of the Jacobian? In other words is the multiplications by the Jacobian determinant built into dolfinx’s building of bi-linear form or should it be added into a potential UFL statement?
The reason why I ask is that I don’t know what dolfinx does exactly and its beyond my bandwidth level to start at this time to go through all of that…
Also its likely to end up a singularity so I’m looking for any comments here on how to resolve that…
