from __future__ import print_function, absolute_import, division #makes KratosMultiphysics backward compatible with python 2.6 and 2.7 # importing the Kratos Library from KratosMultiphysics import * from KratosMultiphysics.ALEApplication import * from KratosMultiphysics.ExternalSolversApplication import * from KratosMultiphysics.IncompressibleFluidApplication import * from KratosMultiphysics.FluidDynamicsApplication import * from KratosMultiphysics.StructuralApplication import * #from KratosMultiphysics.PFEMApplication import PfemUtils from KratosMultiphysics.ULFApplication import * from KratosMultiphysics.MeshingApplication import * import SurfaceTension_monolithic_solver import math import time # Check that KratosMultiphysics was imported in the main script CheckForPreviousImport() # def AddVariables(model_part): # model_part.AddNodalSolutionStepVariable(MESH_DISPLACEMENT) model_part.AddNodalSolutionStepVariable(MESH_VELOCITY) model_part.AddNodalSolutionStepVariable(MESH_REACTION) model_part.AddNodalSolutionStepVariable(MESH_RHS) print("Mesh solver variables added correctly.") # def AddDofs(model_part): for node in model_part.Nodes: node.AddDof(MESH_DISPLACEMENT_X, MESH_REACTION_X) node.AddDof(MESH_DISPLACEMENT_Y, MESH_REACTION_Y) node.AddDof(MESH_DISPLACEMENT_Z, MESH_REACTION_Z) # print("Mesh solver DOFs added correctly.") # import mesh solver base class def CreateSolver(model_part, time_order, reform_dofs_each_step, compute_reactions, box_corner1,box_corner2, domain_size, add_nodes): return MeshSolver(model_part, time_order, reform_dofs_each_step, compute_reactions, compute_reactions, box_corner1,box_corner2, domain_size) class MeshSolver(SurfaceTension_monolithic_solver.STMonolithicSolver): def __init__(self, model_part, time_order, reform_dofs_each_step, compute_reactions, box_corner1,box_corner2, domain_size, add_nodes): #super(MeshSolver, self).__init__(model_part) self.time_order = 2 self.model_part = model_part self.domain_size = 2 self.reform_dofs_each_step = True # neighbour search number_of_avg_elems = 10 number_of_avg_nodes = 10 self.neighbour_search = FindNodalNeighboursProcess(model_part, number_of_avg_elems, number_of_avg_nodes) self.rel_vel_tol = 1E-4 self.abs_vel_tol = 1E-6 self.rel_pres_tol = 1E-4 self.abs_pres_tol= 1E-6 self.echo_level=0 self.alpha = -0.3 self.move_mesh_strategy = 2 self.max_iter = 30 self.divergence_clearance_steps = 0 #self.contact_angle = contact_angle gamma = 1 self.gamma = gamma self.ReformDofSetAtEachStep = True self.CalculateNormDxFlag = False self.MoveMeshFlag = True self.use_slip_conditions = False eul_model_part = 0 self.eul_model_part = eul_model_part self.fluid_neigh_finder = FindNodalNeighboursProcess(model_part,9,18) self.ulf_apply_bc_process = UlfApplyBCProcess(model_part) self.mark_outer_nodes_process = MarkOuterNodesProcess(model_part) bounding_box_corner1_x = -1.00000e+00 bounding_box_corner1_y = -1.00000e+00 bounding_box_corner1_z = -1.00000e+00 bounding_box_corner2_x = 1.01000e+10 bounding_box_corner2_y = 1.01000e+10 bounding_box_corner2_z = 1.01000e+10 box_corner1 = Vector(3); box_corner1[0]=bounding_box_corner1_x; box_corner1[1]=bounding_box_corner1_y; box_corner1[2]=bounding_box_corner1_z; box_corner2 = Vector(3); box_corner2[0]=bounding_box_corner2_x; box_corner2[1]=bounding_box_corner2_y; box_corner2[2]=bounding_box_corner2_z; self.box_corner1 = box_corner1 self.box_corner2 = box_corner2 self.add_nodes=True self.Mesher = TriGenPFEMModeler() self.node_erase_process = NodeEraseProcess(model_part) self.alpha_shape = 3.5 self.condition_neigh_finder = FindConditionsNeighboursProcess(model_part,2, 10) self.mark_fluid_process = MarkFluidProcess(model_part) self.compute_reactions = True self.UlfUtils = UlfUtils() MoveMeshFlag = True # definition of the solvers tol = 1e-8 max_it = 1000 verbosity = 1 m = 10 ##self.linear_solver = AMGCLSolver(AMGCLSmoother.GAUSS_SEIDEL, AMGCLIterativeSolverType.BICGSTAB, tol, max_it, verbosity, m) #pILUPrecond = ILU0Preconditioner() #self.linear_solver = BICGSTABSolver(1e-9, 300, pILUPrecond) #self.linear_solver = SuperLUSolver() # definition of the solvers pILUPrecond = ILU0Preconditioner() self.linear_solver = BICGSTABSolver(1e-9, 300) #def AdaptMesh(self): #pass def Initialize(self): super(MeshSolver, self).Initialize() (self.neighbour_search).Execute() # Initialize mesh model part #for node in self.model_part.Nodes: #zero = Vector(3) #zero[0] = 0.0 #zero[1] = 0.0 #zero[2] = 0.0 #node.SetSolutionStepValue(MESH_REACTION,0,zero) #node.SetSolutionStepValue(MESH_DISPLACEMENT,0,zero) #node.SetSolutionStepValue(MESH_RHS,0,zero) self.solver = StructuralMeshMovingStrategy(self.model_part, self.linear_solver, self.time_order, self.reform_dofs_each_step, self.compute_reactions) #def SetEchoLevel(self, level): (self.solver).SetEchoLevel(0) #def Solve(self): #if(self.reform_dofs_each_step): #(self.neighbour_search).Execute() #(self.solver).Solve() ##(self.solver).MoveMesh() #def MoveMesh(self): #(self.solver).MoveMesh() def GetFluidSolver(self): return super(MeshSolver, self) def GetMeshMotionSolver(self): return self.solver def Solve(self): self.GetMeshMotionSolver().Solve() self.GetFluidSolver().Solve() ## def MoveMesh(self): (self.solver).MoveMesh() #def UpdateReferenceMesh(self): #(self.solver).UpdateReferenceMesh()