diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d646835 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.pyc +__pycache__/ diff --git a/README.md b/README.md index 8885b8b..863753b 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,7 @@ An example mesh showing quadrilateral elements and node numbering: - Python - Parallel Domain Decomposition -- Finite Element concepts +- Finite Element concepts with [xara](https://xara.so) ## File Structure diff --git a/ddm.py b/ddm.py index 5574e27..9768457 100644 --- a/ddm.py +++ b/ddm.py @@ -1,36 +1,35 @@ import numpy as np import pickle from multiprocessing import Process, Queue -import openseespy.opensees as ops -# ----------------------------------------------------------------------------- -# Domain – master process: spawns workers, assembles Schur system, shuts down -# ----------------------------------------------------------------------------- + class Domain: def __init__(self, partitions): """ + master process: spawns workers, assembles Schur system, shuts down + partitions: list of PartitionBuilder instances """ self.partitions = partitions self.workers = [] # list of Process - self.to_subdomain = [] # list of Queue (to worker) - self.from_subdomain = [] # list of Queue (from worker) + self._send = [] # list of Queue (to worker) + self._recv = [] # list of Queue (from worker) self.interface_guess = {} # {(node, dof): value} - self.interface_order = [] # [(node, dof), …] in global order + self.interface_order = [] # [(node, dof), ...] in global order # Spawn one Subdomain worker per partition for part in partitions: - to_q = Queue() - from_q = Queue() + q_send = Queue() + q_recv = Queue() p = Process( target=Domain._subdomain_worker, - args=(pickle.dumps(part), to_q, from_q) + args=(pickle.dumps(part), q_send, q_recv) ) p.start() self.workers.append(p) - self.to_subdomain.append(to_q) - self.from_subdomain.append(from_q) + self._send.append(q_send) + self._recv.append(q_recv) def step(self): """ @@ -43,20 +42,23 @@ def step(self): "dirichlet": self.interface_guess.copy(), "neumann": {} } - for q in self.to_subdomain: + for q in self._send: q.put({"cmd": "solve", "data": interface_data}) - return [q.get() for q in self.from_subdomain] + + return [q.get() for q in self._recv] + def schur_update(self): """ Assemble and solve the global Schur system. """ - for q in self.to_subdomain: + for q in self._send: q.put({"cmd": "schur"}) + S_total = None g_total = None - for q in self.from_subdomain: + for q in self._recv: S_i, g_i = q.get() if S_total is None: S_total = S_i.copy() @@ -76,80 +78,105 @@ def unpack_interface_vector(self, flat_u): } def shutdown(self): - for q in self.to_subdomain: + for q in self._send: q.put("TERMINATE") for p in self.workers: p.join() + @staticmethod - def _subdomain_worker(serialized_partition, to_q, from_q): - import pickle - from ddm import Subdomain - from partitions import LeftQuadPartition, RightQuadPartition + def _subdomain_worker(serialized_partition, q_send, q_recv): partition = pickle.loads(serialized_partition) sd = Subdomain(partition) while True: - msg = to_q.get() + msg = q_send.get() cmd = msg['cmd'] if isinstance(msg, dict) and 'cmd' in msg else msg if cmd == "schur": S_i, g_i = sd.get_schur_data() - from_q.put((S_i, g_i)) + q_recv.put((S_i, g_i)) + elif cmd in ("exit", "TERMINATE"): break + else: raise RuntimeError(f"Unknown message to subdomain worker: {msg}") -# ----------------------------------------------------------------------------- -# Subdomain – local model and Schur contribution -# ----------------------------------------------------------------------------- + class Subdomain: def __init__(self, partition): self.partition = partition + if not hasattr(self.partition, 'node_tags'): raise AttributeError("Partition object lacks node_tags") - ops.model("BasicBuilder", "-ndm", 2, "-ndf", 2) - partition.populate(ops) - ops.system("FullGeneral") - ops.numberer("Plain") - ops.constraints("Plain") - ops.integrator("LoadControl", 1.0) - ops.algorithm("Linear") - ops.analysis("Static") + import opensees.openseespy as _ops + model = _ops.Model("BasicBuilder", "-ndm", 2, "-ndf", 2) + self.model = model + + partition.populate(model) + + model.system("FullGeneral") + model.numberer("Plain") + model.constraints("Plain") + model.integrator("LoadControl", 1.0) + model.algorithm("Linear") + model.analysis("Static") dofs = partition.get_dof_partition() - self.interior_dofs = dofs["interior"] + self.interior_dofs = dofs["interior"] self.interface_dofs = dofs["interface"] + + def apply_interface_conditions(self, interface_data): for (node, dof), val in interface_data.get("dirichlet", {}).items(): - ops.sp(node, dof, val) + self.model.sp(node, dof, val) + for (node, dof), val in interface_data.get("neumann", {}).items(): - ops.load(node, dof, val) + self.model.load(node, dof, val) + def get_schur_data(self): print("Assembling tangent") - rc = ops.analyze(1) + self.model.system("FullGeneral") + self.model.numberer("Plain") + self.model.constraints("Plain") + self.model.test("FixedNumIter", 1) + self.model.integrator("LoadControl", 1.0) + self.model.algorithm("Linear") + self.model.analysis("Static") + + rc = self.model.analyze(1, operation="increment") if rc != 0: raise RuntimeError(f"OpenSees assemble failed (rc={rc})") - print("Assembly complete, fetching K") - N = ops.systemSize() - K = np.array(ops.printA("-ret"), dtype=float).reshape((N, N)) + +# N = self.model.systemSize() + K = self.model.getTangent() # Build global residual - R_global = np.zeros(N) - free_node_list = [n for n in self.partition.node_tags if n not in self.partition.fixed_nodes] - free_nodes = {n: i*2 for i, n in enumerate(free_node_list)} - for node, (Fx, Fy) in self.partition.get_nodal_loads().items(): - if node in free_nodes: - base = free_nodes[node] - R_global[base] = Fx - R_global[base+1] = Fy + R_global = self.model.getResidual() + +# free_node_list = [ +# n for n in self.partition.node_tags +# if n not in self.partition.fixed_nodes +# ] +# free_nodes = {n: i*2 for i, n in enumerate(free_node_list)} + +# for node, (Fx, Fy) in self.partition.get_nodal_loads().items(): +# if node in free_nodes: +# base = free_nodes[node] +# R_global[base] = Fx +# R_global[base+1] = Fy + + +# print(R_global) +# print(f"getResid = {self.model.getResidual()}") I, G = self.interior_dofs, self.interface_dofs + K_II = K[np.ix_(I, I)] K_IG = K[np.ix_(I, G)] K_GI = K[np.ix_(G, I)] @@ -158,7 +185,26 @@ def get_schur_data(self): R_G = R_global[G] KII_inv_KIG = np.linalg.solve(K_II, K_IG) - KII_inv_RI = np.linalg.solve(K_II, R_I) + KII_inv_RI = np.linalg.solve(K_II, R_I) S_local = K_GG - K_GI @ KII_inv_KIG - g_local = R_G - K_GI @ KII_inv_RI + g_local = R_G - K_GI @ KII_inv_RI return S_local, g_local + + +if __name__ == "__main__": + from partitions import LeftQuadPartition, RightQuadPartition + left = LeftQuadPartition() + right = RightQuadPartition() + + print("Starting domain setup") + domain = Domain([left, right]) + domain.interface_order = [(2,1), (2,2), (5,1), (5,2), (8,1), (8,2)] + domain.interface_guess = {key: 0.0 for key in domain.interface_order} + + print("Calling schur_update") + u_gamma = domain.schur_update() + + print("schur_update completed, interface_guess:", domain.interface_guess) + domain.shutdown() + print("Domain shut down") + diff --git a/main.py b/main.py index 7cf829f..2b32f6a 100644 --- a/main.py +++ b/main.py @@ -1,15 +1,18 @@ from ddm import Domain from partitions import LeftQuadPartition, RightQuadPartition -left = LeftQuadPartition() -right = RightQuadPartition() -print("Starting domain setup") -domain = Domain([left, right]) -print("Domain initialized") -domain.interface_order = [(2,1), (2,2), (5,1), (5,2), (8,1), (8,2)] -domain.interface_guess = {key: 0.0 for key in domain.interface_order} -print("Calling schur_update") -u_gamma = domain.schur_update() -print("schur_update completed, interface_guess:", domain.interface_guess) -domain.shutdown() -print("Domain shut down") +if __name__ == "__main__": + left = LeftQuadPartition() + right = RightQuadPartition() + + print("Starting domain setup") + domain = Domain([left, right]) + domain.interface_order = [(2,1), (2,2), (5,1), (5,2), (8,1), (8,2)] + domain.interface_guess = {key: 0.0 for key in domain.interface_order} + + print("Calling schur_update") + u_gamma = domain.schur_update() + + print("schur_update completed, interface_guess:", domain.interface_guess) + domain.shutdown() + print("Domain shut down") diff --git a/matmul b/matmul deleted file mode 100755 index ac29594..0000000 Binary files a/matmul and /dev/null differ diff --git a/model.py b/model.py index 15adecb..c353f7b 100644 --- a/model.py +++ b/model.py @@ -1,4 +1,4 @@ -from openseespy.opensees import * +from opensees.openseespy import * # Define the 2D model with 2 DOFs per node model("BasicBuilder", "-ndm", 2, "-ndf", 2) @@ -56,4 +56,4 @@ # Print the results print("Displacements at interface nodes (node, dof): value") for (node, dof), disp in displacements.items(): - print(f"Node {node}, DOF {dof}: {disp:.12e}") \ No newline at end of file + print(f"Node {node}, DOF {dof}: {disp:.12e}") diff --git a/partitions.py b/partitions.py index afab683..72043d3 100644 --- a/partitions.py +++ b/partitions.py @@ -12,6 +12,7 @@ def __init__(self): self.fixed_nodes = [1] self.nodal_loads = {4: [0.0, -2.0e3]} + def populate(self, model): # 1) Nodes model.node(1, 0.0, 0.0) @@ -32,9 +33,8 @@ def populate(self, model): model.fix(1, 1, 1) # 5) Loads - model.timeSeries("Constant", 1) - model.pattern("Plain", 1, 1) - model.load(4, 0.0, -2.0e3) + model.pattern("Plain", 1, "Constant") + model.load(4, 0.0, -2.0e3, pattern=1) def get_nodal_loads(self): return self.nodal_loads @@ -42,11 +42,15 @@ def get_nodal_loads(self): def get_dof_partition(self): free_nodes = [n for n in self.node_tags if n not in self.fixed_nodes] node_to_dof = {n: i * 2 for i, n in enumerate(free_nodes)} - interior_tags = [t for t in [1, 4, 7] if t not in self.fixed_nodes] + interior_tags = [t for t in [1, 4, 7] if t not in self.fixed_nodes] interface_tags = [t for t in [2, 5, 8] if t not in self.fixed_nodes] interior_dofs = [node_to_dof[t] + d for t in interior_tags for d in (0,1)] interface_dofs = [node_to_dof[t] + d for t in interface_tags for d in (0,1)] - return {"interior": interior_dofs, "interface": interface_dofs} + + return { + "interior": interior_dofs, + "interface": interface_dofs + } # ------------------------------------------------------------------ class RightQuadPartition: