organize, use getResidual to get loads.

.gitignore

@@ -0,0 +1,2 @@
+*.pyc
+__pycache__/

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
 

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")
+

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")

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}")
+    print(f"Node {node}, DOF {dof}: {disp:.12e}")

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,21 +33,24 @@ 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
 
     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: