595 assert(mProblemDimension==DIM+1);
598 unsigned internal_nodes_2d[3] = {3,4,5};
599 unsigned neighbouring_vertices_2d[3][2] = { {1,2}, {2,0}, {0,1} };
602 unsigned internal_nodes_3d[6] = {4,5,6,7,8,9};
603 unsigned neighbouring_vertices_3d[6][2] = { {0,1}, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} };
605 unsigned num_internal_nodes_per_element = DIM==2 ? 3 : 6;
609 iter != mrQuadMesh.GetElementIteratorEnd();
612 for (
unsigned i=0; i<num_internal_nodes_per_element; i++)
614 unsigned global_index;
620 global_index = iter->GetNodeGlobalIndex( internal_nodes_2d[i] );
621 unsigned vertex_0_global_index =iter->GetNodeGlobalIndex( neighbouring_vertices_2d[i][0] );
622 unsigned vertex_1_global_index =iter->GetNodeGlobalIndex( neighbouring_vertices_2d[i][1] );
623 left_val = mCurrentSolution[mProblemDimension*vertex_0_global_index + DIM];
624 right_val = mCurrentSolution[mProblemDimension*vertex_1_global_index + DIM];
628 global_index = iter->GetNodeGlobalIndex( internal_nodes_3d[i] );
629 unsigned vertex_0_global_index =iter->GetNodeGlobalIndex( neighbouring_vertices_3d[i][0] );
630 unsigned vertex_1_global_index =iter->GetNodeGlobalIndex( neighbouring_vertices_3d[i][1] );
631 left_val = mCurrentSolution[mProblemDimension*vertex_0_global_index + DIM];
632 right_val = mCurrentSolution[mProblemDimension*vertex_1_global_index + DIM];
636 mCurrentSolution[mProblemDimension*global_index + DIM] = 0.5 * (left_val + right_val);
674 std::vector<unsigned> rows;
675 std::vector<double> values;
678 bool applySymmetrically = (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY) && symmetricProblem;
680 if (applySymmetrically)
682 if (mDirichletBoundaryConditionsVector ==
nullptr)
684 VecDuplicate(mResidualVector, &mDirichletBoundaryConditionsVector);
695 for (
unsigned i=0; i<mrProblemDefinition.rGetDirichletNodes().size(); i++)
697 unsigned node_index = mrProblemDefinition.rGetDirichletNodes()[i];
699 for (
unsigned j=0; j<DIM; j++)
701 double dirichlet_val = mrProblemDefinition.rGetDirichletNodeValues()[i](j);
706 unsigned dof_index = mProblemDimension*node_index+j;
708 if (type == LINEAR_PROBLEM)
719 val = mCurrentSolution[dof_index] - dirichlet_val;
721 rows.push_back(dof_index);
722 values.push_back(val);
731 if (applySymmetrically)
734 for (
unsigned i=0; i<rows.size(); i++)
736 unsigned col = rows[i];
737 double minus_value = -values[i];
758 if (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY)
763 if (applySymmetrically)
780 if (type!=LINEAR_PROBLEM)
782 for (
unsigned i=0; i<rows.size(); i++)
788 if (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY)
790 for (
unsigned i=0; i<rows.size(); i++)
800 assert(mCompressibilityType==INCOMPRESSIBLE);
803 VecGetOwnershipRange(mResidualVector, &lo, &hi);
805 for (
unsigned i=0; i<mrQuadMesh.GetNumNodes(); i++)
807 if (mrQuadMesh.GetNode(i)->IsInternal())
809 unsigned row = (DIM+1)*i + DIM;
810 if (lo <= (
int)row && (int)row < hi)
812 if (type!=LINEAR_PROBLEM)
816 if (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY)
818 double rhs_vector_val = type==LINEAR_PROBLEM ? 0.0 : mCurrentSolution[row]-0.0;
832 Vec template_vec = mrQuadMesh.GetDistributedVectorFactory()->CreateVec(mProblemDimension);
837 VecDuplicate(template_vec, &mResidualVector);
838 VecDuplicate(mResidualVector, &mLinearSystemRhsVector);
841 mDirichletBoundaryConditionsVector =
nullptr;
849 VecGetOwnershipRange(mResidualVector, &lo, &hi);
856 unsigned num_non_zeros = std::min(75u, mNumDofs);
858 PetscTools::SetupMat(mSystemLhsMatrix, mNumDofs, mNumDofs, num_non_zeros, local_size, local_size);
859 PetscTools::SetupMat(mPreconditionMatrix, mNumDofs, mNumDofs, num_non_zeros, local_size, local_size);
868 int* num_non_zeros_each_row =
new int[mNumDofs];
869 for (
unsigned i=0; i<mNumDofs; i++)
871 num_non_zeros_each_row[i] = 0;
875 iter != mrQuadMesh.GetNodeIteratorEnd();
881 unsigned num_non_zeros_upper_bound = 4 + 30*iter->GetNumContainingElements();
883 num_non_zeros_upper_bound = std::min(num_non_zeros_upper_bound, mNumDofs);
885 unsigned i = iter->GetIndex();
887 num_non_zeros_each_row[mProblemDimension*i + 0] = num_non_zeros_upper_bound;
888 num_non_zeros_each_row[mProblemDimension*i + 1] = num_non_zeros_upper_bound;
889 num_non_zeros_each_row[mProblemDimension*i + 2] = num_non_zeros_upper_bound;
891 if (mCompressibilityType==INCOMPRESSIBLE)
893 if (!iter->IsInternal())
895 num_non_zeros_each_row[mProblemDimension*i + 3] = num_non_zeros_upper_bound;
899 num_non_zeros_each_row[mProblemDimension*i + 3] = 1;
918#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2)
919 MatCreate(PETSC_COMM_WORLD,local_size,local_size,mNumDofs,mNumDofs,&mSystemLhsMatrix);
920 MatCreate(PETSC_COMM_WORLD,local_size,local_size,mNumDofs,mNumDofs,&mPreconditionMatrix);
922 MatCreate(PETSC_COMM_WORLD,&mSystemLhsMatrix);
923 MatCreate(PETSC_COMM_WORLD,&mPreconditionMatrix);
924 MatSetSizes(mSystemLhsMatrix,local_size,local_size,mNumDofs,mNumDofs);
925 MatSetSizes(mPreconditionMatrix,local_size,local_size,mNumDofs,mNumDofs);
930 MatSetType(mSystemLhsMatrix, MATSEQAIJ);
931 MatSetType(mPreconditionMatrix, MATSEQAIJ);
932 MatSeqAIJSetPreallocation(mSystemLhsMatrix, PETSC_DEFAULT, num_non_zeros_each_row);
933 MatSeqAIJSetPreallocation(mPreconditionMatrix, PETSC_DEFAULT, num_non_zeros_each_row);
937 int* num_non_zeros_each_row_in_diag =
new int[local_size];
938 int* num_non_zeros_each_row_off_diag =
new int[local_size];
939 for (
unsigned i=0; i<
unsigned(local_size); i++)
941 num_non_zeros_each_row_in_diag[i] = num_non_zeros_each_row[lo+i];
942 num_non_zeros_each_row_off_diag[i] = num_non_zeros_each_row[lo+i];
944 if (num_non_zeros_each_row_in_diag[i] > local_size)
946 num_non_zeros_each_row_in_diag[i] = local_size;
950 MatSetType(mSystemLhsMatrix, MATMPIAIJ);
951 MatSetType(mPreconditionMatrix, MATMPIAIJ);
952 MatMPIAIJSetPreallocation(mSystemLhsMatrix, PETSC_DEFAULT, num_non_zeros_each_row_in_diag, PETSC_DEFAULT, num_non_zeros_each_row_off_diag);
953 MatMPIAIJSetPreallocation(mPreconditionMatrix, PETSC_DEFAULT, num_non_zeros_each_row_in_diag, PETSC_DEFAULT, num_non_zeros_each_row_off_diag);
956 MatSetFromOptions(mSystemLhsMatrix);
957 MatSetFromOptions(mPreconditionMatrix);
958#if (PETSC_VERSION_MAJOR == 3)
959 MatSetOption(mSystemLhsMatrix, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
960 MatSetOption(mPreconditionMatrix, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
962 MatSetOption(mSystemLhsMatrix, MAT_IGNORE_OFF_PROC_ENTRIES);
963 MatSetOption(mPreconditionMatrix, MAT_IGNORE_OFF_PROC_ENTRIES);
973 delete [] num_non_zeros_each_row;