Merge pull request #1982 from Hukai0/draft/kyle

WIP(cgraph):add linear_elasticity_eigen_3d_example

Author: weihuayi

.gitignore

@@ -75,3 +75,4 @@ tags
 .DS_Store
 *env
 tmp
+*.msh

docs/Gemfile.lock

@@ -85,7 +85,7 @@ GEM
     rb-fsevent (0.11.0)
     rb-inotify (0.10.1)
       ffi (~> 1.0)
-    rexml (3.3.9)
+    rexml (3.4.2)
     rouge (3.26.0)
     safe_yaml (1.0.5)
     sassc (2.4.0)

example/cfd/bam_break.py

@@ -0,0 +1,37 @@
+from fealpy.backend import backend_manager as bm
+import matplotlib.pyplot as plt
+bm.set_backend("numpy")
+from fealpy.mesh.node_mesh import BamBreak
+from fealpy.cfd.simulation.sph.particle_solver_new import ParticleSystem, BamBreakSolver, Visualizer
+
+dx = 0.03
+dy = 0.03
+maxstep = 2000
+
+# 创建网格
+mesh = BamBreak.from_bam_break_domain(dx, dy)
+
+# 创建粒子系统
+particle_system = ParticleSystem(dx, dy)
+particle_system.initialize_particles(mesh.node, mesh.nodedata)
+
+# 创建SPH求解器
+sph_solver = BamBreakSolver(particle_system)
+
+# 初始可视化
+color = bm.where(particle_system.particles["tag"], "red", "blue")
+plt.scatter(particle_system.particles["position"][:, 0], 
+            particle_system.particles["position"][:, 1], 
+            c=color, s=5)
+plt.grid(True)
+plt.show()
+
+# 运行模拟
+sph_solver.run_simulation(maxstep, draw_interval=10, reinitalize_interval=30)
+
+# 最终结果可视化
+visualizer = Visualizer(particle_system)
+visualizer.final_plot()
+
+
+

example/cfd/test/stationary_stokes.py

@@ -52,7 +52,7 @@
     help = "Type of refinement strategy, default is uniform_refine")
 
 parser.add_argument('--maxit',
-    default = 5, type = int,
+    default = 3, type = int,
     help = "Maximum number of iterations for the solver, default is 5")
 
 parser.add_argument('--maxstep',

example/cgraph/dipole_example.py

@@ -0,0 +1,49 @@
+
+import json
+import fealpy.cgraph as cgraph
+
+WORLD_GRAPH = cgraph.WORLD_GRAPH
+
+pde = cgraph.create("DipoleAntenna3D")
+mesher = cgraph.create("Dipole3d")
+spacer = cgraph.create("FunctionSpace")
+uher = cgraph.create("FEFunction")
+dipole_antenna_eq = cgraph.create("DipoleAntennaEquation")
+dbc = cgraph.create("DirichletBC")
+solver = cgraph.create("IterativeSolver")
+
+
+# 连接节点
+mesher(cr = 0.05, sr0 = 1.9, sr1 = 2.4, L = 2.01, G = 0.01)
+spacer(type = "first_nedelec", mesh=mesher().mesh, p=1)
+dipole_antenna_eq(
+    space=spacer(),
+    q=3,
+    diffusion = pde().diffusion,
+    reaction = pde().reaction,
+    source = pde().source,
+    Y = pde().Y,
+    ID = mesher().ID1
+)
+
+dbc(
+    gd=pde().dirichlet,
+    isDDof=mesher().ID2,
+    form=dipole_antenna_eq().operator,
+    F=dipole_antenna_eq().source
+)
+
+solver(
+    A=dbc().A,
+    b=dbc().F,
+    solver = "minres"
+)
+
+
+# 最终连接到图输出节点上
+WORLD_GRAPH.output(mesh=mesher().mesh, uh=solver())
+
+WORLD_GRAPH.register_error_hook(print)
+WORLD_GRAPH.execute()
+print(WORLD_GRAPH.get())
+

example/cgraph/dld_microfluidic_chip_2d_example.py

@@ -0,0 +1,39 @@
+import json
+import fealpy.cgraph as cgraph
+
+
+WORLD_GRAPH = cgraph.WORLD_GRAPH
+
+pde = cgraph.create("DLDMicrofluidicChip2D")
+mesher = cgraph.create("DLDMicrofluidicChipMesh2d")
+uspacer = cgraph.create("TensorFunctionSpace")
+pspacer = cgraph.create("FunctionSpace")
+dld_eq = cgraph.create("StokesEquation")
+solver = cgraph.create("CGSolver")
+postprocess = cgraph.create("VPDecoupling")
+
+mesher(lc = 0.02)
+uspacer(mesh = mesher(), p=2, gd = 2)
+pspacer(mesh = mesher(), p=1)
+pde(radius = mesher().radius,
+    centers = mesher().centers,
+    inlet_boundary = mesher().inlet_boundary,
+    outlet_boundary = mesher().outlet_boundary,
+    wall_boundary = mesher().wall_boundary)
+dld_eq(uspace = uspacer(),
+    pspace = pspacer(),
+    velocity_dirichlet = pde().velocity_dirichlet, 
+    pressure_dirichlet = pde().pressure_dirichlet, 
+    is_velocity_boundary = pde().is_velocity_boundary, 
+    is_pressure_boundary = pde().is_pressure_boundary)
+solver(A = dld_eq().bform,
+       b = dld_eq().lform)
+postprocess(out = solver().out, 
+            uspace = uspacer(),
+            mesh = mesher())
+
+# 最终连接到图输出节点上
+WORLD_GRAPH.output(uh = postprocess().uh,u_x = postprocess().u_x, u_y = postprocess().u_y, ph = postprocess().ph)
+WORLD_GRAPH.error_listeners.append(print)
+WORLD_GRAPH.execute()
+print(WORLD_GRAPH.get())

example/cgraph/dld_microfluidic_chip_3d_example.py

@@ -0,0 +1,40 @@
+import json
+import fealpy.cgraph as cgraph
+
+
+WORLD_GRAPH = cgraph.WORLD_GRAPH
+
+pde = cgraph.create("DLDMicrofluidicChip3D")
+mesher = cgraph.create("DLDMicrofluidicChipMesh3d")
+uspacer = cgraph.create("TensorFunctionSpace")
+pspacer = cgraph.create("FunctionSpace")
+dld_eq = cgraph.create("StokesEquation")
+solver = cgraph.create("CGSolver")
+postprocess = cgraph.create("VPDecoupling")
+
+mesher(lc = 0.04)
+uspacer(mesh = mesher(), p=2, gd = 3)
+pspacer(mesh = mesher(), p=1)
+pde(thickness = mesher().thickness,
+    radius = mesher().radius,
+    centers = mesher().centers,
+    inlet_boundary = mesher().inlet_boundary,
+    outlet_boundary = mesher().outlet_boundary,
+    wall_boundary = mesher().wall_boundary)
+dld_eq(uspace = uspacer(),
+    pspace = pspacer(),
+    velocity_dirichlet = pde().velocity_dirichlet, 
+    pressure_dirichlet = pde().pressure_dirichlet, 
+    is_velocity_boundary = pde().is_velocity_boundary, 
+    is_pressure_boundary = pde().is_pressure_boundary)
+solver(A = dld_eq().bform,
+       b = dld_eq().lform)
+postprocess(out = solver().out, 
+            uspace = uspacer(),
+            mesh = mesher())
+
+# 最终连接到图输出节点上
+WORLD_GRAPH.output(uh = postprocess().uh)
+WORLD_GRAPH.error_listeners.append(print)
+WORLD_GRAPH.execute()
+print(WORLD_GRAPH.get())

example/cgraph/euler_bernoulli_beam_2d_example.py

@@ -0,0 +1,49 @@
+import fealpy.cgraph as cgraph
+from fealpy.backend import backend_manager as bm
+
+WORLD_GRAPH = cgraph.WORLD_GRAPH
+
+model = cgraph.create("Beam2d")
+mesher = cgraph.create("CreateMesh")
+beam_materialer = cgraph.create("BeamMaterial")
+spacer = cgraph.create("FunctionSpace")
+beam_model = cgraph.create("Beam")
+dbc = cgraph.create("StructuralDirichletBC")
+solver = cgraph.create("DirectSolver")
+postprocess = cgraph.create("UDecoupling")
+
+# 连接节点
+node = bm.array([[0], [5],[7.5]], dtype=bm.float64)
+cell = bm.array([[0, 1],[1,2]] , dtype=bm.int32)
+
+mesher(node = node, cell = cell)
+spacer(type="lagrange", mesh=mesher(), p=1)
+
+beam_materialer(property="Steel", beam_type="Euler-Bernoulli beam", beam_E=200e9, beam_nu=0.3, I=118.6e-6)
+
+beam_model(
+    space = spacer(),
+    beam_E = beam_materialer().E,
+    beam_nu = beam_materialer().nu,
+    I = beam_materialer().I,
+    distributed_load = model().f,
+    beam_type = "euler_bernoulli_2d",
+)
+
+dbc(
+    gd=model().dirichlet,
+    isDDof=model().dirichlet_dof_index,
+    K = beam_model().K,
+    F=beam_model().F,
+    space=spacer()
+)
+
+solver(A = dbc().K,
+       b = dbc().F)
+
+
+# 最终连接到图输出节点上
+WORLD_GRAPH.output(mesh=mesher(), u=solver())
+WORLD_GRAPH.register_error_hook(print)
+WORLD_GRAPH.execute()
+print(WORLD_GRAPH.get())

example/cgraph/helmholtz_2d_example.py

@@ -12,6 +12,7 @@
 helmholtz_eq = cgraph.create("HelmholtzEquation")
 dbc = cgraph.create("DirichletBC")
 solver = cgraph.create("CGSolver")
+mesh2d3d = cgraph.create("MeshDimensionUpgrading")
 
 
 # 连接节点
@@ -35,7 +36,19 @@
     x0=dbc().uh
 )
 # 最终连接到图输出节点上
-WORLD_GRAPH.output(mesh=mesher(), uh=solver())
+WORLD_GRAPH.output(
+    mesh=mesh2d3d(mesh=mesher(), z=solver()),
+    uh=solver()
+)
 
 WORLD_GRAPH.execute()
-print(WORLD_GRAPH.get())
+result = WORLD_GRAPH.get()
+mesh = result["mesh"]
+
+from matplotlib import pyplot as plt
+
+fig = plt.figure()
+axes = fig.add_subplot(projection="3d")
+mesh.add_plot(axes)
+
+plt.show()

example/cgraph/incompressible_cylinder.py

@@ -0,0 +1,73 @@
+import json
+import fealpy.cgraph as cgraph
+
+
+WORLD_GRAPH = cgraph.WORLD_GRAPH
+
+pde = cgraph.create("IncompressibleCylinder2d")
+uspacer = cgraph.create("TensorFunctionSpace")
+pspacer = cgraph.create("FunctionSpace")
+dbc_u = cgraph.create("ProjectDBC")
+dbc_p = cgraph.create("ProjectDBC")
+simulation = cgraph.create("IncompressibleNSIPCS")
+timeline = cgraph.create("CFDTimeline")
+IncompressibleNSRun = cgraph.create("IncompressibleNSIPCSRun")
+
+pde(
+    mu = 0.001,
+    rho = 1.0,
+    cx = 0.5,
+    cy = 0.2,
+    radius = 0.07,
+    n_circle = 100,
+    h = 0.06)
+uspacer(mesh = pde().mesh, p=2, gd = 2)
+pspacer(mesh = pde().mesh, p=1)
+dbc_u(
+    space = uspacer(), 
+    dirichlet = pde().velocity_dirichlet, 
+    is_boundary = pde().is_velocity_boundary
+)
+dbc_p(
+    space = pspacer(), 
+    dirichlet = pde().pressure_dirichlet, 
+    is_boundary = pde().is_pressure_boundary
+)
+simulation(
+    constitutive = 1,
+    mu = pde().mu,
+    rho = pde().rho,
+    source = pde().source,
+    uspace = uspacer(),
+    pspace = pspacer(),
+    is_pressure_boundary = pde().is_pressure_boundary,
+    apply_bcu = dbc_u().apply_bc,
+    apply_bcp = dbc_p().apply_bc,
+    q = 2
+) 
+timeline(
+    T0 = 0.0,
+    T1 = 1.0,
+    NT = 1000
+)
+IncompressibleNSRun(
+    T0=timeline().T0,
+    T1=timeline().T1,
+    NL=timeline().NL,
+    uspace = uspacer(), 
+    pspace = pspacer(), 
+    velocity_0 = pde().velocity_0,
+    pressure_0 = pde().pressure_0,
+    is_pressure_boundary = pde().is_pressure_boundary,
+    predict_velocity = simulation().predict_velocity,
+    correct_pressure = simulation().correct_pressure,
+    correct_velocity = simulation().correct_velocity,
+    mesh = pde().mesh,
+    output_dir = "/home/libz/cylinder"
+)
+
+WORLD_GRAPH.output(uh_x = IncompressibleNSRun().uh_x)
+WORLD_GRAPH.error_listeners.append(print)
+WORLD_GRAPH.execute()
+print(WORLD_GRAPH.get())
+