Merge pull request #13912 from KratosMultiphysics/PR_corrections_laplacian_and_fluids

[IgaApplication] small corrections and typos in Fluid and Laplacian workflow of IgaApplication

Author: andrewgorgi

applications/IgaApplication/custom_conditions/sbm_fluid_condition_dirichlet.cpp

@@ -39,9 +39,9 @@ void SbmFluidConditionDirichlet::CalculateAll(
     KRATOS_TRY
 
     const auto& r_geometry = GetGeometry();
-    const SizeType number_of_nodes = r_geometry.size();
+    const std::size_t number_of_nodes = r_geometry.size();
 
-    const SizeType mat_size = number_of_nodes * (mDim+1);
+    const std::size_t mat_size = number_of_nodes * (mDim+1);
     //resizing as needed the LHS
     if(rLeftHandSideMatrix.size1() != mat_size)
         rLeftHandSideMatrix.resize(mat_size,mat_size,false);
@@ -113,16 +113,25 @@ void SbmFluidConditionDirichlet::CalculateAll(
     stress_old(1, 0) = r_stress_vector[2];      stress_old(1, 1) = r_stress_vector[1];         
     Vector traction_old_iteration = prod(stress_old, n_tensor); // This results in a 2x1 vector.
 
+    Matrix DB_contribution_w = ZeroMatrix(2, 2);
+    Matrix DB_contribution = ZeroMatrix(2, 2);
+
     for (IndexType i = 0; i < number_of_nodes; i++) {
-        for (IndexType j = 0; j < number_of_nodes; j++) {
-            for (IndexType idim = 0; idim < 2; idim++) {
+        for (IndexType idim = 0; idim < 2; idim++) {
+            DB_contribution_w(0, 0) = DB_voigt(0, 2*i+idim);
+            DB_contribution_w(0, 1) = DB_voigt(2, 2*i+idim);
+            DB_contribution_w(1, 0) = DB_voigt(2, 2*i+idim);
+            DB_contribution_w(1, 1) = DB_voigt(1, 2*i+idim);
+
+            for (IndexType j = 0; j < number_of_nodes; j++) {
+                // Compute the traction vector: sigma * n.
+                Vector traction_nitsche_w = prod(DB_contribution_w, n_tensor);
 
                 // Penalty term
                 rLeftHandSideMatrix(3*i+idim, 3*j+idim) += mHsum(0,i)*mHsum(0,j)* penalty_integration;
 
                 for (IndexType jdim = 0; jdim < 2; jdim++) {
                     // Extract the 2x2 block for the control point i from the sigma matrix.
-                    Matrix DB_contribution = ZeroMatrix(2, 2);
                     DB_contribution(0, 0) = DB_voigt(0, 2*j+jdim);
                     DB_contribution(0, 1) = DB_voigt(2, 2*j+jdim);
                     DB_contribution(1, 0) = DB_voigt(2, 2*j+jdim);
@@ -134,12 +143,16 @@ void SbmFluidConditionDirichlet::CalculateAll(
                     rLeftHandSideMatrix(3*i+idim, 3*j+jdim) -= H(0, i) * traction(idim) * integration_weight;
 
                     // skew-symmetric Nitsche term
-                    rLeftHandSideMatrix(3*j+jdim, 3*i+idim) += mHsum(0, i) * traction(idim) * integration_weight;
+                    rLeftHandSideMatrix(3*i+idim, 3*j+jdim) += mHsum(0, j) * traction_nitsche_w(jdim) * integration_weight;
                 }
 
                 // integration by parts PRESSURE
                 rLeftHandSideMatrix(3*i+idim, 3*j+2) += H(0,j)* ( H(0,i) * mNormalParameterSpace[idim] )
                         * integration_weight;
+                
+                // Nitsche term --> q term
+                rLeftHandSideMatrix(3*j+mDim, 3*i+idim) -= H(0,j)* ( mHsum(0,i) * mNormalParameterSpace[idim] )
+                        * integration_weight;
             }
         }
 
@@ -153,14 +166,18 @@ void SbmFluidConditionDirichlet::CalculateAll(
             rRightHandSideVector(3*i+idim) -= pressure_old_iteration * ( H(0,i) * mNormalParameterSpace[idim] ) * integration_weight;
 
             // skew-symmetric Nitsche term
-            Matrix DB_contribution = ZeroMatrix(2, 2); // Extract the 2x2 block for the control point i from the DB_voigt.
             DB_contribution(0, 0) = DB_voigt(0, 2*i+idim); 
             DB_contribution(0, 1) = DB_voigt(2, 2*i+idim); 
             DB_contribution(1, 0) = DB_voigt(2, 2*i+idim);
             DB_contribution(1, 1) = DB_voigt(1, 2*i+idim); 
             // Compute the traction vector: sigma * n.
             Vector traction = prod(DB_contribution, n_tensor);
-            rRightHandSideVector(3*i+idim) -= velocity_old_iteration[idim] * traction(idim) * integration_weight;
+            for (IndexType jdim = 0; jdim < mDim; jdim++) {
+                rRightHandSideVector(3*i+idim) -= velocity_old_iteration[jdim] * traction(jdim) * integration_weight;
+            }
+            // Nitsche term --> q term
+            rRightHandSideVector(3*i+mDim) += velocity_old_iteration[idim] * ( H(0,i) * mNormalParameterSpace[idim] )
+                        * integration_weight;
         }
     }
 
@@ -184,7 +201,11 @@ void SbmFluidConditionDirichlet::CalculateAll(
             // Compute the traction vector: sigma * n.
             Vector traction = prod(sigma_block, n_tensor);
             // skew-symmetric Nitsche term
-            rRightHandSideVector[3*i+idim] += u_D[idim] * traction(idim) * integration_weight;
+            for (IndexType jdim = 0; jdim < mDim; jdim++) {
+                rRightHandSideVector[3*i+idim] += u_D[jdim] * traction(jdim) * integration_weight;
+            }
+            // Nitsche term --> q term
+            rRightHandSideVector[3*i+mDim] -= u_D[idim] * H(0,i)*mNormalParameterSpace[idim] * integration_weight;
 
         }
     }
@@ -230,7 +251,7 @@ void SbmFluidConditionDirichlet::InitializeMemberVariables()
 void SbmFluidConditionDirichlet::InitializeSbmMemberVariables()
 {
     const auto& r_geometry = this->GetGeometry();
-    const SizeType number_of_nodes = r_geometry.size();
+    const std::size_t number_of_nodes = r_geometry.size();
 
     // Retrieve projection
     Condition candidate_closest_skin_segment_1 = this->GetValue(NEIGHBOUR_CONDITIONS)[0] ;
@@ -285,8 +306,8 @@ void SbmFluidConditionDirichlet::CalculateB(
         Matrix& rB, 
         const ShapeDerivativesType& r_DN_DX) const
 {
-    const SizeType number_of_control_points = GetGeometry().size();
-    const SizeType mat_size = number_of_control_points * 2; // Only 2 DOFs per node in 2D
+    const std::size_t number_of_control_points = GetGeometry().size();
+    const std::size_t mat_size = number_of_control_points * 2; // Only 2 DOFs per node in 2D
 
     // Resize B matrix to 3 rows (strain vector size) and appropriate number of columns
     if (rB.size1() != 3 || rB.size2() != mat_size)
@@ -311,7 +332,7 @@ void SbmFluidConditionDirichlet::ApplyConstitutiveLaw(
         ConstitutiveLaw::Parameters& rValues,
         ConstitutiveVariables& rConstitutiveVariables) const
 {
-    const SizeType number_of_nodes = GetGeometry().size();
+    const std::size_t number_of_nodes = GetGeometry().size();
 
     // Set constitutive law flags:
     Flags& ConstitutiveLawOptions=rValues.GetOptions();
@@ -355,7 +376,7 @@ int SbmFluidConditionDirichlet::Check(const ProcessInfo& rCurrentProcessInfo) co
 void SbmFluidConditionDirichlet::EquationIdVector(EquationIdVectorType &rResult, const ProcessInfo &rCurrentProcessInfo) const
 {
     const GeometryType& rGeom = this->GetGeometry();
-    const SizeType number_of_control_points = GetGeometry().size();
+    const std::size_t number_of_control_points = GetGeometry().size();
     const unsigned int LocalSize = (mDim + 1) * number_of_control_points;
 
     if (rResult.size() != LocalSize)
@@ -379,7 +400,7 @@ void SbmFluidConditionDirichlet::GetDofList(
 {
     KRATOS_TRY;
 
-    const SizeType number_of_control_points = GetGeometry().size();
+    const std::size_t number_of_control_points = GetGeometry().size();
 
     rElementalDofList.resize(0);
     rElementalDofList.reserve((mDim + 1) * number_of_control_points);
@@ -397,8 +418,8 @@ void SbmFluidConditionDirichlet::GetDofList(
 void SbmFluidConditionDirichlet::GetSolutionCoefficientVector(
         Vector& rValues) const
 {
-    const SizeType number_of_control_points = GetGeometry().size();
-    const SizeType mat_size = number_of_control_points * mDim;
+    const std::size_t number_of_control_points = GetGeometry().size();
+    const std::size_t mat_size = number_of_control_points * mDim;
 
     if (rValues.size() != mat_size)
         rValues.resize(mat_size, false);

applications/IgaApplication/custom_conditions/sbm_fluid_condition_dirichlet.h

@@ -35,7 +35,6 @@ namespace Kratos
         KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(SbmFluidConditionDirichlet);
 
         /// Size types
-        typedef std::size_t SizeType;
         typedef std::size_t IndexType;
 
         /// Type for shape function derivatives container

applications/IgaApplication/custom_conditions/sbm_laplacian_condition_neumann.cpp

@@ -39,71 +39,20 @@ void SbmLaplacianConditionNeumann::CalculateLocalSystem(
 {
     KRATOS_TRY
 
-    ConvectionDiffusionSettings::Pointer p_settings = rCurrentProcessInfo[CONVECTION_DIFFUSION_SETTINGS];
-    const auto& r_unknown_var = p_settings->GetUnknownVariable();
+    const std::size_t mat_size = GetGeometry().size() * 1;
 
-    const auto& r_geometry = this->GetGeometry();
-    const SizeType number_of_nodes = r_geometry.PointsNumber();
-    if (rRightHandSideVector.size() != number_of_nodes) {
-        rRightHandSideVector.resize(number_of_nodes, false);
-    }
-    if (rLeftHandSideMatrix.size1() != number_of_nodes || rLeftHandSideMatrix.size2() != number_of_nodes) {
-        rLeftHandSideMatrix.resize(number_of_nodes, number_of_nodes, false);
-    }
-
-    noalias(rRightHandSideVector) = ZeroVector(number_of_nodes);
-    noalias(rLeftHandSideMatrix) = ZeroMatrix(number_of_nodes, number_of_nodes);
-
-    // Integration
-    const auto& r_integration_points = r_geometry.IntegrationPoints();
-    const auto& r_DN_De = r_geometry.ShapeFunctionsLocalGradients(r_geometry.GetDefaultIntegrationMethod());
-    
-    // Initialize DN_DX
-    Matrix DN_DX(number_of_nodes,mDim);
+    if (rRightHandSideVector.size() != mat_size)
+        rRightHandSideVector.resize(mat_size);
+    noalias(rRightHandSideVector) = ZeroVector(mat_size);
 
-    // Calculating the PHYSICAL SPACE derivatives (it is avoided storing them to minimize storage)
-    noalias(DN_DX) = r_DN_De[0]; // prod(r_DN_De[0],InvJ0);
-    
-    const Matrix& H = r_geometry.ShapeFunctionsValues();
+    if (rLeftHandSideMatrix.size1() != mat_size)
+        rLeftHandSideMatrix.resize(mat_size, mat_size);
+    noalias(rLeftHandSideMatrix) = ZeroMatrix(mat_size, mat_size);
 
-    // compute gradient Taylor expansion contribution: grad_H_sum
-    Matrix grad_H_sum = ZeroMatrix(3, number_of_nodes);
-    ComputeGradientTaylorExpansionContribution(grad_H_sum);
-    
-    // dot product grad cdot n
-    Matrix H_grad_N_dot_n = ZeroMatrix(1, number_of_nodes);
-    noalias(row(H_grad_N_dot_n, 0)) = row(grad_H_sum, 0) * mTrueNormal[0] + row(grad_H_sum, 1) * mTrueNormal[1] + row(grad_H_sum, 2) * mTrueNormal[2];
+    CalculateLeftHandSide(rLeftHandSideMatrix,rCurrentProcessInfo);
 
-    Matrix DN_dot_n_tilde = ZeroMatrix(1, number_of_nodes);
-    for (IndexType i = 0; i < number_of_nodes; ++i)
-    {
-        // grad N cdot n_tilde
-        for (IndexType idim = 0; idim < mDim; idim++) {
-            DN_dot_n_tilde(0, i)  += DN_DX(i, idim) * mNormalParameterSpace[idim];
-        } 
-    }
-    
-    // Assembly
-    // compute Neumann contributions
-    noalias(rLeftHandSideMatrix) += prod(trans(H), H_grad_N_dot_n) * mTrueDotSurrogateNormal * r_integration_points[0].Weight(); // * std::abs(determinant_jacobian_vector[point_number]) ;
-    noalias(rLeftHandSideMatrix) -= prod(trans(H), DN_dot_n_tilde)                           * r_integration_points[0].Weight() ; // * std::abs(DetJ0) ;
+    CalculateRightHandSide(rRightHandSideVector,rCurrentProcessInfo);
 
-    Vector t_N(number_of_nodes);
-    for (IndexType i = 0; i < number_of_nodes; ++i)
-    {
-        t_N[i] = mpProjectionNode->GetValue(HEAT_FLUX);
-    }
-    // Neumann Contributions
-    noalias(rRightHandSideVector) += prod(prod(trans(H), H), t_N) * mTrueDotSurrogateNormal * r_integration_points[0].Weight(); // * std::abs(determinant_jacobian_vector[point_number]);
-
-    Vector temp(number_of_nodes);
-    // RHS = ExtForces - K*temp;
-    for (IndexType i = 0; i < number_of_nodes; i++) {
-        temp[i] = r_geometry[i].GetSolutionStepValue(r_unknown_var);
-    }
-    // RHS -= K*temp
-    noalias(rRightHandSideVector) -= prod(rLeftHandSideMatrix,temp);
-
     KRATOS_CATCH("")
 
 }
@@ -198,7 +147,7 @@ void SbmLaplacianConditionNeumann::CalculateLeftHandSide(
 {
     ConvectionDiffusionSettings::Pointer p_settings = rCurrentProcessInfo[CONVECTION_DIFFUSION_SETTINGS];
     const auto& r_geometry = this->GetGeometry();
-    const SizeType number_of_nodes = r_geometry.PointsNumber();
+    const std::size_t number_of_nodes = r_geometry.PointsNumber();
 
     if (rLeftHandSideMatrix.size1() != number_of_nodes || rLeftHandSideMatrix.size2() != number_of_nodes) {
         rLeftHandSideMatrix.resize(number_of_nodes, number_of_nodes, false);
@@ -246,7 +195,9 @@ void SbmLaplacianConditionNeumann::CalculateRightHandSide(
     const ProcessInfo& rCurrentProcessInfo)
 {
     const auto& r_geometry = this->GetGeometry();
-    const SizeType number_of_nodes = r_geometry.PointsNumber();
+    ConvectionDiffusionSettings::Pointer p_settings = rCurrentProcessInfo[CONVECTION_DIFFUSION_SETTINGS];
+    const auto& r_unknown_var = p_settings->GetUnknownVariable();
+    const std::size_t number_of_nodes = r_geometry.PointsNumber();
     if (rRightHandSideVector.size() != number_of_nodes) {
         rRightHandSideVector.resize(number_of_nodes, false);
     }
@@ -265,12 +216,24 @@ void SbmLaplacianConditionNeumann::CalculateRightHandSide(
     }
     // Neumann Contributions
     noalias(rRightHandSideVector) += prod(prod(trans(H), H), t_N) * mTrueDotSurrogateNormal * r_integration_points[0].Weight(); // * std::abs(determinant_jacobian_vector[point_number]);
+
+    // --- Corresponding RHS ---
+    Matrix left_hand_side = ZeroMatrix(number_of_nodes);
+    CalculateLeftHandSide(left_hand_side,rCurrentProcessInfo);
+    Vector temp(number_of_nodes);
+    // RHS = ExtForces - K*temp;
+    for (IndexType i = 0; i < number_of_nodes; i++) {
+        temp[i] = r_geometry[i].GetSolutionStepValue(r_unknown_var);
+    }
+    // RHS -= K*temp
+    noalias(rRightHandSideVector) -= prod(left_hand_side,temp);
+
 }
 
 void SbmLaplacianConditionNeumann::ComputeGradientTaylorExpansionContribution(Matrix& grad_H_sum)
 {
     const auto& r_geometry = this->GetGeometry();
-    const SizeType number_of_nodes = r_geometry.PointsNumber();
+    const std::size_t number_of_nodes = r_geometry.PointsNumber();
     const auto& r_DN_De = r_geometry.ShapeFunctionsLocalGradients(r_geometry.GetDefaultIntegrationMethod());
 
     // Compute all the derivatives of the basis functions involved
@@ -384,7 +347,7 @@ void SbmLaplacianConditionNeumann::EquationIdVector(
     const auto& r_unknown_var = p_settings->GetUnknownVariable();
 
     const auto& r_geometry = GetGeometry();
-    const SizeType number_of_nodes = r_geometry.size();
+    const std::size_t number_of_nodes = r_geometry.size();
 
     if (rResult.size() !=  number_of_nodes)
         rResult.resize(number_of_nodes, false);
@@ -403,7 +366,7 @@ void SbmLaplacianConditionNeumann::GetDofList(
     const auto& r_unknown_var = p_settings->GetUnknownVariable();
 
     const auto& r_geometry = GetGeometry();
-    const SizeType number_of_nodes = r_geometry.size();
+    const std::size_t number_of_nodes = r_geometry.size();
 
     rElementalDofList.resize(0);
     rElementalDofList.reserve(number_of_nodes);