36 #ifndef ABSTRACTCONTINUUMMECHANICSSOLVER_HPP_
37 #define ABSTRACTCONTINUUMMECHANICSSOLVER_HPP_
39 #include "ContinuumMechanicsProblemDefinition.hpp"
40 #include "CompressibilityType.hpp"
41 #include "LinearSystem.hpp"
42 #include "OutputFileHandler.hpp"
43 #include "ReplicatableVector.hpp"
44 #include "AbstractTetrahedralMesh.hpp"
45 #include "QuadraticMesh.hpp"
46 #include "DistributedQuadraticMesh.hpp"
47 #include "Warnings.hpp"
48 #include "PetscException.hpp"
49 #include "GaussianQuadratureRule.hpp"
51 #include "MechanicsEventHandler.hpp"
52 #include "CommandLineArguments.hpp"
53 #include "VtkMeshWriter.hpp"
61 typedef enum _ApplyDirichletBcsType
64 NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY,
65 NONLINEAR_PROBLEM_APPLY_TO_EVERYTHING
66 } ApplyDirichletBcsType;
86 template<
unsigned DIM>
300 std::string outputDirectory,
301 CompressibilityType compressibilityType);
347 void CreateVtkOutput(std::string spatialSolutionName=
"Spatial solution");
375 template<
unsigned DIM>
378 std::string outputDirectory,
379 CompressibilityType compressibilityType)
380 : mrQuadMesh(rQuadMesh),
381 mrProblemDefinition(rProblemDefinition),
382 mOutputDirectory(outputDirectory),
383 mpOutputFileHandler(nullptr),
384 mpQuadratureRule(nullptr),
385 mpBoundaryQuadratureRule(nullptr),
386 mCompressibilityType(compressibilityType),
387 mResidualVector(nullptr),
388 mSystemLhsMatrix(nullptr),
389 mPreconditionMatrix(nullptr)
391 assert(DIM==2 || DIM==3);
397 if ((p_quad_mesh ==
nullptr) && (p_distributed_quad_mesh ==
nullptr))
399 EXCEPTION(
"Continuum mechanics solvers require a quadratic mesh");
429 template<
unsigned DIM>
432 if (mpOutputFileHandler)
434 delete mpOutputFileHandler;
437 if (mpQuadratureRule)
439 delete mpQuadratureRule;
440 delete mpBoundaryQuadratureRule;
451 if (mDirichletBoundaryConditionsVector)
457 template<
unsigned DIM>
459 std::string fileExtension,
470 std::stringstream file_name;
472 file_name << fileName;
473 if (counterToAppend >= 0)
475 file_name <<
"_" << counterToAppend;
477 file_name <<
"." << fileExtension;
479 out_stream p_file = mpOutputFileHandler->OpenOutputFile(file_name.str());
481 std::vector<c_vector<double,DIM> >& r_spatial_solution = rGetSpatialSolution();
482 for (
unsigned i=0; i<r_spatial_solution.size(); i++)
488 for (
unsigned j=0; j<DIM; j++)
490 *p_file << r_spatial_solution[i](j) <<
" ";
499 template<
unsigned DIM>
508 std::stringstream file_name;
510 file_name <<
"pressure";
511 if (counterToAppend >= 0)
513 file_name <<
"_" << counterToAppend;
517 out_stream p_file = mpOutputFileHandler->OpenOutputFile(file_name.str());
519 std::vector<double> &r_pressure = rGetPressures();
520 for (
unsigned i = 0; i < r_pressure.size(); i++)
522 for (
unsigned j = 0; j < DIM; j++)
524 *p_file << mrQuadMesh.GetNode(i)->rGetLocation()[j] <<
" ";
527 *p_file << r_pressure[i] <<
"\n";
535 template<
unsigned DIM>
538 if (writeOutput && (mOutputDirectory==
""))
540 EXCEPTION(
"Can't write output if no output directory was given in constructor");
542 mWriteOutput = writeOutput;
546 template<
unsigned DIM>
549 if (this->mOutputDirectory==
"")
551 EXCEPTION(
"No output directory was given so no output was written, cannot convert to VTK format");
556 mesh_writer.
AddPointData(spatialSolutionName, this->rGetSpatialSolution());
558 if (mCompressibilityType==INCOMPRESSIBLE)
564 std::vector<double> element_attribute;
566 iter != this->mrQuadMesh.GetElementIteratorEnd();
569 element_attribute.push_back(iter->GetAttribute());
571 mesh_writer.
AddCellData(
"Attribute", element_attribute);
577 template<
unsigned DIM>
580 assert(mProblemDimension==DIM+1);
582 mPressureSolution.clear();
583 mPressureSolution.resize(mrQuadMesh.GetNumNodes());
585 for (
unsigned i=0; i<mrQuadMesh.GetNumNodes(); i++)
587 mPressureSolution[i] = mCurrentSolution[mProblemDimension*i + DIM];
589 return mPressureSolution;
592 template<
unsigned DIM>
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);
671 template<
unsigned DIM>
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)
778 if (type!=LINEAR_PROBLEM)
780 for (
unsigned i=0; i<rows.size(); i++)
786 if (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY)
788 for (
unsigned i=0; i<rows.size(); i++)
795 template<
unsigned DIM>
798 assert(mCompressibilityType==INCOMPRESSIBLE);
801 VecGetOwnershipRange(mResidualVector, &lo, &hi);
803 for (
unsigned i=0; i<mrQuadMesh.GetNumNodes(); i++)
805 if (mrQuadMesh.GetNode(i)->IsInternal())
807 unsigned row = (DIM+1)*i + DIM;
808 if (lo <= (
int)row && (int)row < hi)
810 if (type!=LINEAR_PROBLEM)
814 if (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY)
816 double rhs_vector_val = type==LINEAR_PROBLEM ? 0.0 : mCurrentSolution[row]-0.0;
827 template<
unsigned DIM>
830 Vec template_vec = mrQuadMesh.GetDistributedVectorFactory()->CreateVec(mProblemDimension);
835 VecDuplicate(template_vec, &mResidualVector);
836 VecDuplicate(mResidualVector, &mLinearSystemRhsVector);
839 mDirichletBoundaryConditionsVector =
nullptr;
847 VecGetOwnershipRange(mResidualVector, &lo, &hi);
854 unsigned num_non_zeros = std::min(75u, mNumDofs);
856 PetscTools::SetupMat(mSystemLhsMatrix, mNumDofs, mNumDofs, num_non_zeros, local_size, local_size);
857 PetscTools::SetupMat(mPreconditionMatrix, mNumDofs, mNumDofs, num_non_zeros, local_size, local_size);
866 int* num_non_zeros_each_row =
new int[mNumDofs];
867 for (
unsigned i=0; i<mNumDofs; i++)
869 num_non_zeros_each_row[i] = 0;
873 iter != mrQuadMesh.GetNodeIteratorEnd();
879 unsigned num_non_zeros_upper_bound = 4 + 30*iter->GetNumContainingElements();
881 num_non_zeros_upper_bound = std::min(num_non_zeros_upper_bound, mNumDofs);
883 unsigned i = iter->GetIndex();
885 num_non_zeros_each_row[mProblemDimension*i + 0] = num_non_zeros_upper_bound;
886 num_non_zeros_each_row[mProblemDimension*i + 1] = num_non_zeros_upper_bound;
887 num_non_zeros_each_row[mProblemDimension*i + 2] = num_non_zeros_upper_bound;
889 if (mCompressibilityType==INCOMPRESSIBLE)
891 if (!iter->IsInternal())
893 num_non_zeros_each_row[mProblemDimension*i + 3] = num_non_zeros_upper_bound;
897 num_non_zeros_each_row[mProblemDimension*i + 3] = 1;
916 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
917 MatCreate(PETSC_COMM_WORLD,local_size,local_size,mNumDofs,mNumDofs,&mSystemLhsMatrix);
918 MatCreate(PETSC_COMM_WORLD,local_size,local_size,mNumDofs,mNumDofs,&mPreconditionMatrix);
920 MatCreate(PETSC_COMM_WORLD,&mSystemLhsMatrix);
921 MatCreate(PETSC_COMM_WORLD,&mPreconditionMatrix);
922 MatSetSizes(mSystemLhsMatrix,local_size,local_size,mNumDofs,mNumDofs);
923 MatSetSizes(mPreconditionMatrix,local_size,local_size,mNumDofs,mNumDofs);
928 MatSetType(mSystemLhsMatrix, MATSEQAIJ);
929 MatSetType(mPreconditionMatrix, MATSEQAIJ);
930 MatSeqAIJSetPreallocation(mSystemLhsMatrix, PETSC_DEFAULT, num_non_zeros_each_row);
931 MatSeqAIJSetPreallocation(mPreconditionMatrix, PETSC_DEFAULT, num_non_zeros_each_row);
935 int* num_non_zeros_each_row_in_diag =
new int[local_size];
936 int* num_non_zeros_each_row_off_diag =
new int[local_size];
937 for (
unsigned i=0; i<
unsigned(local_size); i++)
939 num_non_zeros_each_row_in_diag[i] = num_non_zeros_each_row[lo+i];
940 num_non_zeros_each_row_off_diag[i] = num_non_zeros_each_row[lo+i];
942 if (num_non_zeros_each_row_in_diag[i] > local_size)
944 num_non_zeros_each_row_in_diag[i] = local_size;
948 MatSetType(mSystemLhsMatrix, MATMPIAIJ);
949 MatSetType(mPreconditionMatrix, MATMPIAIJ);
950 MatMPIAIJSetPreallocation(mSystemLhsMatrix, PETSC_DEFAULT, num_non_zeros_each_row_in_diag, PETSC_DEFAULT, num_non_zeros_each_row_off_diag);
951 MatMPIAIJSetPreallocation(mPreconditionMatrix, PETSC_DEFAULT, num_non_zeros_each_row_in_diag, PETSC_DEFAULT, num_non_zeros_each_row_off_diag);
954 MatSetFromOptions(mSystemLhsMatrix);
955 MatSetFromOptions(mPreconditionMatrix);
956 #if (PETSC_VERSION_MAJOR == 3) //PETSc 3.x.x
957 MatSetOption(mSystemLhsMatrix, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
958 MatSetOption(mPreconditionMatrix, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
960 MatSetOption(mSystemLhsMatrix, MAT_IGNORE_OFF_PROC_ENTRIES);
961 MatSetOption(mPreconditionMatrix, MAT_IGNORE_OFF_PROC_ENTRIES);
971 delete [] num_non_zeros_each_row;
974 #endif // ABSTRACTCONTINUUMMECHANICSSOLVER_HPP_
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
void ApplyDirichletBoundaryConditions(ApplyDirichletBcsType type, bool symmetricProblem)
GaussianQuadratureRule< DIM > * mpQuadratureRule
void WriteCurrentPressureSolution(int counterToAppend=-1)
Vec mLinearSystemRhsVector
AbstractTetrahedralMesh< DIM, DIM > & mrQuadMesh
virtual std::vector< c_vector< double, DIM > > & rGetSpatialSolution()=0
void AllocateMatrixMemory()
void WriteCurrentSpatialSolution(std::string fileName, std::string fileExtension, int counterToAppend=-1)
#define EXCEPTION(message)
virtual ~AbstractContinuumMechanicsSolver()
std::vector< double > & rGetPressures()
Vec mDirichletBoundaryConditionsVector
void AddCellData(std::string name, std::vector< double > data)
std::vector< c_vector< double, DIM > > mSpatialSolution
virtual unsigned GetNumNodes() const
AbstractContinuumMechanicsSolver(AbstractTetrahedralMesh< DIM, DIM > &rQuadMesh, ContinuumMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory, CompressibilityType compressibilityType)
void AddIdentityBlockForDummyPressureVariables(ApplyDirichletBcsType type)
bool OptionExists(const std::string &rOption)
void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)
void RemovePressureDummyValuesThroughLinearInterpolation()
std::vector< double > & rGetCurrentSolution()
ContinuumMechanicsProblemDefinition< DIM > & mrProblemDefinition
NodeIterator GetNodeIteratorBegin(bool skipDeletedNodes=true)
OutputFileHandler * mpOutputFileHandler
void CreateVtkOutput(std::string spatialSolutionName="Spatial solution")
CompressibilityType mCompressibilityType
std::vector< double > mCurrentSolution
std::vector< double > mPressureSolution
unsigned mProblemDimension
static CommandLineArguments * Instance()
GaussianQuadratureRule< DIM-1 > * mpBoundaryQuadratureRule
std::string mOutputDirectory
void AddPointData(std::string name, std::vector< double > data)
void SetWriteOutput(bool writeOutput=true)