Add Pfem Melting application

applications/ConvectionDiffusionApplication/custom_elements/eulerian_conv_diff.cpp

@@ -112,6 +112,7 @@ namespace Kratos
             }
         }
 
+
         //Some auxilary definitions
         BoundedMatrix<double,TNumNodes, TNumNodes> aux1 = ZeroMatrix(TNumNodes, TNumNodes); //terms multiplying dphi/dt
         BoundedMatrix<double,TNumNodes, TNumNodes> aux2 = ZeroMatrix(TNumNodes, TNumNodes); //terms multiplying phi
@@ -129,20 +130,21 @@ namespace Kratos
                  for(unsigned int k=0; k<TDim; k++)
                     vel_gauss[k] += N[i]*(Variables.v[i][k]*Variables.theta + Variables.vold[i][k]*(1.0-Variables.theta));
             }
+
             const double norm_vel = norm_2(vel_gauss);
             array_1d<double, TNumNodes > a_dot_grad = prod(DN_DX, vel_gauss);
 
             const double tau = this->CalculateTau(Variables,norm_vel,h);
 
             //terms multiplying dphi/dt (aux1)
-            noalias(aux1) += (1.0+tau*Variables.beta*Variables.div_v)*outer_prod(N, N);
+            noalias(aux1) += (1.0+tau*Variables.beta*Variables.div_v)* 0.25 * IdentityMatrix(4, 4); //outer_prod(N, N); //0.25 * IdentityMatrix(4, 4);
             noalias(aux1) +=  tau*outer_prod(a_dot_grad, N);
 
             //terms which multiply the gradient of phi
             noalias(aux2) += (1.0+tau*Variables.beta*Variables.div_v)*outer_prod(N, a_dot_grad);
             noalias(aux2) += tau*outer_prod(a_dot_grad, a_dot_grad);
-        }
 
+        }
         //adding the second and third term in the formulation
         noalias(rLeftHandSideMatrix)  = (Variables.dt_inv*Variables.density*Variables.specific_heat + Variables.theta*Variables.beta*Variables.div_v)*aux1;
         noalias(rRightHandSideVector) = (Variables.dt_inv*Variables.density*Variables.specific_heat - (1.0-Variables.theta)*Variables.beta*Variables.div_v)*prod(aux1,Variables.phi_old);
@@ -261,13 +263,17 @@ namespace Kratos
 				  rVariables.v[i] = GetGeometry()[i].FastGetSolutionStepValue(rVelocityVar);
 				  rVariables.vold[i] = GetGeometry()[i].FastGetSolutionStepValue(rVelocityVar,1);
 				  //active_convection=true;
+
+
+
 			}
 
 			if (IsDefinedMeshVelocityVariable)
             {
 				  const Variable<array_1d<double, 3 > >& rMeshVelocityVar = my_settings->GetMeshVelocityVariable();
 				  rVariables.v[i] -= GetGeometry()[i].FastGetSolutionStepValue(rMeshVelocityVar);
 				  rVariables.vold[i] -= GetGeometry()[i].FastGetSolutionStepValue(rMeshVelocityVar,1);
+
 				  //active_convection=true;
 			}
 

applications/FluidDynamicsApplication/custom_elements/vms.h

@@ -100,7 +100,7 @@ namespace Kratos
  */
 template< unsigned int TDim,
           unsigned int TNumNodes = TDim + 1 >
-class VMS : public Element
+class KRATOS_API(FLUID_DYNAMICS_APPLICATION) VMS : public Element
 {
 public:
     ///@name Type Definitions

applications/FluidDynamicsApplication/custom_strategies/schemes/residualbased_predictorcorrector_velocity_bossak_scheme_turbulent.h

@@ -329,7 +329,7 @@ namespace Kratos {
                             array_1d<double, 3 > & OldVelocity = (itNode)->FastGetSolutionStepValue(VELOCITY, 1);
 
                             noalias(itNode->FastGetSolutionStepValue(MESH_VELOCITY)) = itNode->FastGetSolutionStepValue(VELOCITY);
-                            UpdateDisplacement(CurrentDisplacement, OldDisplacement, OldVelocity, OldAcceleration, CurrentAcceleration);
+                            //UpdateDisplacement(CurrentDisplacement, OldDisplacement, OldVelocity, OldAcceleration, CurrentAcceleration);
                         }
                         else
                         {
@@ -407,7 +407,7 @@ namespace Kratos {
                   if((itNode)->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) < 1e-15)
 			{
 			    noalias(itNode->FastGetSolutionStepValue(MESH_VELOCITY)) = itNode->FastGetSolutionStepValue(VELOCITY);
-			    UpdateDisplacement(CurrentDisplacement, OldDisplacement, OldVelocity, OldAcceleration, CurrentAcceleration);
+			    //UpdateDisplacement(CurrentDisplacement, OldDisplacement, OldVelocity, OldAcceleration, CurrentAcceleration);
 			}
 			else
 			{

applications/MeshingApplication/CMakeLists.txt

@@ -18,6 +18,7 @@ if(${USE_TETGEN_NONFREE_TPL} MATCHES ON)
     include_directories(${TETGEN_INCLUDE})
     add_library(tetgen STATIC ${TETGEN_INCLUDE}/tetgen.cxx ${TETGEN_INCLUDE}/predicates.cxx)
     target_compile_definitions(tetgen PRIVATE -DTETLIBRARY)
+    add_definitions(-DUSE_TETGEN_NONFREE_TPL)
 endif(${USE_TETGEN_NONFREE_TPL} MATCHES ON)
 
 ## generate variables with the sources

applications/MeshingApplication/custom_python/add_meshers_to_python.cpp

@@ -28,7 +28,7 @@
     #include "external_includes/tetgen_cdt.h"
 #else
     #define REAL double
-#endif 
+#endif
 
 #if USE_TRIANGLE_NONFREE_TPL
     #include "external_includes/trigen_pfem_refine.h"
@@ -37,7 +37,7 @@
     #include "external_includes/trigen_glass_forming.h"
     #include "external_includes/trigen_droplet_refine.h"
     #include "external_includes/trigen_cdt.h"
-#endif 
+#endif
 
 namespace Kratos
 {
@@ -200,9 +200,11 @@ void TriRegenerateMeshVMS(TriGenPFEMModelerVMS& Mesher, char* ElementName, char*
 
 void  AddMeshersToPython(pybind11::module& m)
 {
-
+
+#pragma message "A ****************************************************************************************"
 #ifdef USE_TETGEN_NONFREE_TPL
     // Class that allows 3D adaptive remeshing (inserting and erasing nodes)
+    #pragma message "B ****************************************************************************************"
     py::class_<TetGenPfemModeler, TetGenPfemModeler::Pointer >(m, "TetGenPfemModeler")
     .def(py::init< >())
     .def("ReGenerateMesh",TetRegenerateMesh)
@@ -218,7 +220,7 @@ void  AddMeshersToPython(pybind11::module& m)
     .def(py::init< >())
     .def("ReGenerateMesh",TetRegenerateMeshContact)
     ;
-    
+
     py::class_<TetGenCDT, TetGenCDT::Pointer >(m, "TetGenCDT")
     .def(py::init< >())
     .def("GenerateCDT",GenerateCDT)
@@ -229,7 +231,7 @@ void  AddMeshersToPython(pybind11::module& m)
     .def("ReGenerateMesh",&TetGenPfemModelerVms::ReGenerateMesh)
     ;
 #endif
-    
+
 #if USE_TRIANGLE_NONFREE_TPL
     // Class that allows 2D adaptive remeshing (inserting and erasing nodes)
     py::class_<TriGenPFEMModeler, TriGenPFEMModeler::Pointer >(m, "TriGenPFEMModeler")

applications/MeshingApplication/external_includes/tetgen_pfem_refine.h

@@ -156,7 +156,7 @@ namespace Kratos
 			PointVector res(max_results);
 			DistanceVector res_distances(max_results);
 			Node<3> work_point(0,0.0,0.0,0.0);
-			KRATOS_WATCH(h_factor)
+			//KRATOS_WATCH(h_factor)
  			//if the remove_node switch is activated, we check if the nodes got too close
 			if (rem_nodes==true)
 			{
@@ -238,7 +238,7 @@ namespace Kratos
 
 
 			double first_part_time = auxiliary.elapsed();
-			std::cout << "mesh generation time = " << first_part_time << std::endl;
+			//std::cout << "mesh generation time = " << first_part_time << std::endl;
 
 			//generate Kratos Tetrahedra3D4
 			int el_number = out.numberoftetrahedra;
@@ -406,7 +406,7 @@ namespace Kratos
 				}
 
 			}
-			std::cout << "time for passing alpha shape" << alpha_shape_time.elapsed() << std::endl;
+			//std::cout << "time for passing alpha shape" << alpha_shape_time.elapsed() << std::endl;
 
 			//freeing unnecessary memory
 			in.deinitialize();
@@ -426,7 +426,7 @@ namespace Kratos
 				in2.pointlist[base+1] = (nodes_begin + i)->Y();
 				in2.pointlist[base+2] = (nodes_begin + i)->Z();
 			}
-			std::cout << "qui" << std::endl;
+			//std::cout << "qui" << std::endl;
 			in2.numberoftetrahedra = number_of_preserved_elems;
 			in2.tetrahedronlist = new int[in2.numberoftetrahedra * 4];
 			in2.tetrahedronvolumelist = new double[in2.numberoftetrahedra];
@@ -516,7 +516,7 @@ namespace Kratos
 
 			//q - creates quality mesh, with the default radius-edge ratio set to 2.0
 
- 			std::cout << "mesh recreation time" << mesh_recreation_time.elapsed() << std::endl;
+ 			//std::cout << "mesh recreation time" << mesh_recreation_time.elapsed() << std::endl;
 
 
 			//PAVEL
@@ -577,7 +577,7 @@ namespace Kratos
 				}
 			}
 
-			std::cout << "During refinement we added " << outnew.numberofpoints-n_points_before_refinement<< "nodes " <<std::endl;
+			//std::cout << "During refinement we added " << outnew.numberofpoints-n_points_before_refinement<< "nodes " <<std::endl;
 
 
 			bucket_size = 50;
@@ -758,7 +758,7 @@ ModelPart::NodesContainerType& ModelNodes = ThisModelPart.Nodes();
 				(ThisModelPart.Elements()).push_back(p_element);
 
 			}
-			std::cout << "time for adding elems" << adding_elems.elapsed() << std::endl;;
+			//std::cout << "time for adding elems" << adding_elems.elapsed() << std::endl;;
 			ThisModelPart.Elements().Sort();
 
 			boost::timer adding_neighb;
@@ -782,7 +782,7 @@ ModelPart::NodesContainerType& ModelNodes = ThisModelPart.Nodes();
 						neighb(i) = Element::WeakPointer();
 				}
 			}
-			std::cout << "time for adding neigbours" << adding_neighb.elapsed() << std::endl;;
+			//std::cout << "time for adding neigbours" << adding_neighb.elapsed() << std::endl;;
 
 
 
@@ -797,7 +797,7 @@ ModelPart::NodesContainerType& ModelNodes = ThisModelPart.Nodes();
 				in->FastGetSolutionStepValue(IS_BOUNDARY) = 0;
 			}
 
-                        std::cout << "reset the boundary flag" << adding_neighb.elapsed() << std::endl;;
+                        //std::cout << "reset the boundary flag" << adding_neighb.elapsed() << std::endl;;
 			//***********************************************************************************
 			//***********************************************************************************
 			boost::timer adding_faces;
@@ -844,7 +844,7 @@ ModelPart::NodesContainerType& ModelNodes = ThisModelPart.Nodes();
 			}
 			outnew.deinitialize();
 			outnew.initialize();
-			std::cout << "time for adding faces" << adding_faces.elapsed() << std::endl;;
+			//std::cout << "time for adding faces" << adding_faces.elapsed() << std::endl;;
 
 
 
@@ -870,7 +870,7 @@ ModelPart::NodesContainerType& ModelNodes = ThisModelPart.Nodes();
 
 
 			double second_part_time = auxiliary.elapsed();
-			std::cout << "second part time = " << second_part_time - first_part_time << std::endl;
+			//std::cout << "second part time = " << second_part_time - first_part_time << std::endl;
 
 
 			KRATOS_CATCH("")