How to use the Cutting Application (Serial and Parallel)

From KratosWiki
Jump to: navigation, search

In many 3D models that are calculated, the part of the domain that we're interested in is a small fraction of the full domain and therefore in the postprocess we create cuts of this small part, leaving all the rest unused. This is ok as long as the domain is small, but when it has a million or more elements, we end up with GBs of data that makes reading results a slow task.

x x

The idea behind this utility is to print a postprocess file that contains only the results of cutting planes that have been defined as input. Therefore results in the full domain are discarded once they are not (internally) needed anymore by the selected solver scheme. User can choose the layer in which each plane will be printed, to turning them on and off is trivial.


To do this we will have to add the next tasks in our python file:

  • Create a new, empty model part that will later contain the cutting planes and add the same variables to this new model part
  • Define the cutting planes that we want with a point a normal unit vector and call the constructor of the utility so that the nodes and triangles are created.
  • at each time step that we want to print results, call the function that interpolates the results from the original tetrahedra into the triangles.
  • Tell Kratos to print the results here instead of the other model part.

Contents

Serial(only one processor)

1. Importing the Meshing Application and defining the new model part

since this utility is part of the Meshing App, we'll have to import it. For example in a fluid problem:

       ...
       #importing applications
       import applications_interface
       applications_interface.Import_IncompressibleFluidApplication = True
       applications_interface.Import_FluidDynamicsApplication = True
       applications_interface.Import_MeshingApplication = True 		        #import the meshing app interface
       applications_interface.ImportApplications(kernel, kratos_applications_path)
       
       ## from now on the order is not anymore crucial
       ##################################################################
       ##################################################################
       from KratosIncompressibleFluidApplication import *
       from KratosFluidDynamicsApplication import *
       from KratosMeshingApplication import * 						#import the meshing app
       
       fluid_model_part = ModelPart("FluidPart");                                      #original model part   
       cut_model_part = ModelPart("CutPart");						#define a new model part in which the results will be printed
       ...


In the default scripts (for example, the ones generated by GID) , only elements are printed. We will have to change this because the mesh we are about to generate contains conditions instead of elements: (if you get an empty mesh as output it might be that your cutting planes are not cutting the domain or that you forgot to change this)

       write_conditions = WriteConditionsFlag.WriteConditions 			# ORIGINAL: write_conditions = WriteConditionsFlag.WriteElementsOnly
       gid_io = GidIO(input_file_name,gid_mode,multifile,deformed_mesh_flag, write_conditions)

2. Defining the cutting planes and calling the utility

Both the points and normal vectors that define the planes are standard vectors:

       normal_X = Vector([1.0  ,  0.0 ,   0.0])
       normal_Y = Vector([0.0  ,  1.0 ,   0.0])
       normal_Z = Vector([0.0  ,  0.0 ,   1.0])
       
       X_0c0  = Vector([0.0  ,  0.0 ,   0.0])
       X_0c5  = Vector([0.5  ,  0.0 ,   0.0])
       X_1c0  = Vector([1.0  ,  0.0 ,   0.0])
       X_1c5  = Vector([1.5  ,  0.0 ,   0.0])
       X_2c0  = Vector([2.0  ,  0.0 ,   0.0])
       X_2c5  = Vector([2.5  ,  0.0 ,   0.0])
       X_3c0  = Vector([3.0  ,  0.0 ,   0.0])
       Z_0c0  = Vector([0.0005  ,  0.0005 ,   0.0005])
       Y_0c0  = Vector([0.0005  ,  0.0005 ,   0.0005])

We call the constructor:

       #constructor of the cutting application (this is done only once)
       Cut_App = Cutting_Application();

And now we find the smallest edge. This is not necessary but recommended to set the tolerance (explained below). If this is not done, absolute tolerance is assumed. If this value is small enough there will be no problems. However, if it happens to be larger than hald the smallest edge you're likely to have problems and Kratos might crash.

       #find the smallest edge (can be done only once if  the mesh is fixed)
       Cut_App.FindSmallestEdge(fluid_model_part)

We are ready to fill the new model part with the cutting planes. If skin conditions (using triangles) have been used, for example an inlet, we can also print the in the same way.

       #Generate the cut. This function has to be called once per cutting plane.
       #factor is a double : 0.5<factor<0.0 (if FindSmallestEdge was used) that
       #will be used to decide 
       #whether to keep the position of a mesh node when the cutting point is
       #closer than (tolerance* smallest edge). a high value will generate
       #bigger triangles but the result will be far from a plane. on the other
       #hand a value close to zero will generate lots of small triangles.
       #layer(integer) is the layer in which the generated plane will be printed
       #Cut_App.GenerateCut(fluid_model_part,cut_model_part,versor,Xp, layer , factor)
       
       Cut_App.GenerateCut(fluid_model_part,cut_model_part,normal_X,X_0c0, 1 , 0.01)
       Cut_App.AddSkinConditions(fluid_model_part,cut_model_part, 2)
       Cut_App.GenerateCut(fluid_model_part,cut_model_part,normal_X,X_0c5, 3 , 0.1)
       Cut_App.GenerateCut(fluid_model_part,cut_model_part,normal_X,X_1c0, 4 , 0.1)
       Cut_App.GenerateCut(fluid_model_part,cut_model_part,normal_X,X_1c5, 5 , 0.1)
       Cut_App.GenerateCut(fluid_model_part,cut_model_part,normal_X,X_2c0, 6 , 0.1)
       Cut_App.GenerateCut(fluid_model_part,cut_model_part,normal_X,X_2c5, 7 , 0.1)
       Cut_App.GenerateCut(fluid_model_part,cut_model_part,normal_X,X_3c0, 8 , 0.1)
       Cut_App.GenerateCut(fluid_model_part,cut_model_part,normal_Y,Y_0c0, 9 , 0.01)
       Cut_App.GenerateCut(fluid_model_part,cut_model_part,normal_Z,Z_0c0, 10 , 0.1)

Now that we have the new mesh, we must print it:

       #mesh to be printed
       mesh_name = 0.0
       gid_io.InitializeMesh( mesh_name)
       gid_io.WriteMesh( cut_model_part.GetMesh() ) 					#originally: gid_io.WriteMesh( fluid_model_part.GetMesh() )
       gid_io.FinalizeMesh()
       
       gid_io.InitializeResults(mesh_name,(cut_model_part).GetMesh())			#originally: gid_io.InitializeResults(mesh_name,(fluid_model_part).GetMesh())

3. Updating and priting results

At each time step: To begin with we must clone the time steps:

       fluid_model_part.CloneTimeStep(time)
       cut_model_part.CloneTimeStep(time)

And after the system has been solved, in the timesteps that we want our results to be printed we will have to first interpolate the results. For example if we wanted to print these two variables:

           #gid_io.WriteNodalResults(PRESSURE,fluid_model_part.Nodes,time,0)
           #gid_io.WriteNodalResults(VELOCITY,fluid_model_part.Nodes,time,0)

We call UpdateCutData, that interpolates the results in all the variables in the current timestep, and then we write using the new model part:

           Cut_App.UpdateCutData(cut_model_part,fluid_model_part)
           gid_io.WriteNodalResults(PRESSURE,cut_model_part.Nodes,time,0)
           gid_io.WriteNodalResults(VELOCITY,cut_model_part.Nodes,time,0)

And that's it, the results are printed in the planes.

Parallel

It performs the same task but for domains that are distributed in different processors. The user interface is almost the same. The main difference is that it is included in the Trilinos Aplication, so it will have to be imported instead of the Meshing App:

       ...
       #importing applications
       import applications_interface
       applications_interface.Import_IncompressibleFluidApplication = True
       applications_interface.Import_KratosTrilinosApplication = True
       applications_interface.Import_KratosMetisApplication = True
       applications_interface.ImportApplications(kernel, kratos_applications_path)
       
       ## from now on the order is not anymore crucial
       ##################################################################
       ##################################################################
       from KratosIncompressibleFluidApplication import *
       from KratosTrilinosApplication import *
       from KratosMetisApplication import *
       
       #defining a model part for the fluid and one for the structure
       fluid_model_part = ModelPart("FluidPart")                         #original model part   
       cut_model_part = ModelPart("FluidPart")                           #define a new model part in which the results will be printed
       ...

Other difference is the constructor, which requires a communicator:

       Comm = CreateCommunicator()
       Cut_App = TrilinosCuttingApplication(Comm);


Finally, when printing the results we must have in mind that each of the processes have both its local nodes and the ones belonging to its neighbour. So if we printed the results in all the nodes we would sometimes have duplicated results, which would cause problems when visualizing all the parts toghether in the postprocess. To fix this, we'll write the results only in the LOCAL nodes:

          Cut_App.UpdateCutData(cut_model_part,fluid_model_part)
          local_nodes = cut_model_part.GetCommunicator().LocalMesh().Nodes       #these are only the local nodes
          gid_io.WriteNodalResults(PRESSURE,local_nodes,time,0)   #and so we save the results only in them.
          gid_io.WriteNodalResults(VELOCITY,local_nodes,time,0)

The rest of the steps remain unchanged

Personal tools
Categories