diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index bb9f98e..25efed2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -37,4 +37,4 @@ jobs: - name: "Run for ${{ matrix.os }}" shell: bash working-directory: ${{runner.workspace}} - run: ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" -D CTEST_CONFIGURATION_TYPE=RelWithDebInfo -D LABELS_FOR_SUBPROJECTS="gsElasticity" -D CTEST_SITE="${{ matrix.os }}_[actions]" -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON" -D GISMO_OPTIONAL="${{ github.event.repository.name }}" -Q + run: ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" -D CTEST_CONFIGURATION_TYPE=RelWithDebInfo -D LABELS_FOR_SUBPROJECTS="gsElasticity-examples" -D CTEST_SITE="${{ matrix.os }}_[actions]" -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON" -D GISMO_OPTIONAL="${{ github.event.repository.name }}" -Q diff --git a/CMakeLists.txt b/CMakeLists.txt index 3be35e8..372bf81 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -54,7 +54,6 @@ add_definitions(-DELAST_DATA_DIR="${CMAKE_CURRENT_SOURCE_DIR}/filedata/") # add example files if(GISMO_BUILD_EXAMPLES) - add_custom_target(${PROJECT_NAME}) add_custom_target(${PROJECT_NAME}-examples) add_subdirectory(examples) else() diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index db004e2..033026b 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -11,16 +11,15 @@ aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR} FILES) foreach(file ${FILES}) add_gismo_executable(${file}) get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}-examples") if(GISMO_BUILD_EXAMPLES) - set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") + set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") else(GISMO_BUILD_EXAMPLES) set_target_properties(${tarname} PROPERTIES - FOLDER "${PROJECT_NAME}" - EXCLUDE_FROM_ALL TRUE) + FOLDER "${PROJECT_NAME}" + EXCLUDE_FROM_ALL TRUE) endif(GISMO_BUILD_EXAMPLES) add_dependencies(${PROJECT_NAME}-examples ${tarname}) - add_dependencies(${PROJECT_NAME} ${tarname}) # install the example executables (optionally) install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) endforeach(file ${FILES}) \ No newline at end of file diff --git a/examples/biceps_activeMuscle_3D.cpp b/examples/biceps_activeMuscle_3D.cpp index d47c098..cda1d0f 100644 --- a/examples/biceps_activeMuscle_3D.cpp +++ b/examples/biceps_activeMuscle_3D.cpp @@ -39,12 +39,12 @@ int main(int argc, char* argv[]){ gsVector<> fiberDirection(3); fiberDirection << 1.,0.,0.; // space discretization - index_t numUniRefDirX = 2; + index_t numUniRefDirX = 0; index_t numUniRef = 0; index_t numDegElev = 0; bool subgridOrTaylorHood = false; // output - index_t numPlotPoints = 10000; + index_t numPlotPoints = 0; // minimalistic user interface for terminal gsCmdLine cmd("This is a simulation of active muscle behavior."); diff --git a/examples/biceps_activeMuscle_3Dt.cpp b/examples/biceps_activeMuscle_3Dt.cpp index 5c4b6fd..78236d8 100644 --- a/examples/biceps_activeMuscle_3Dt.cpp +++ b/examples/biceps_activeMuscle_3Dt.cpp @@ -39,7 +39,7 @@ int main(int argc, char* argv[]){ gsVector<> fiberDirection(3); fiberDirection << 1.,0.,0.; // space discretization - index_t numUniRefDirX = 2; + index_t numUniRefDirX = 0; index_t numUniRef = 0; index_t numDegElev = 0; bool subgridOrTaylorHood = false; @@ -47,7 +47,7 @@ int main(int argc, char* argv[]){ real_t timeSpan = 2; real_t timeStep = 0.1; // output - index_t numPlotPoints = 100000; + index_t numPlotPoints = 0; // minimalistic user interface for terminal gsCmdLine cmd("This is a simulation of active muscle behavior."); diff --git a/examples/flappingBeam_FSI2_coupledFSI_2Dt.cpp b/examples/flappingBeam_FSI2_coupledFSI_2Dt.cpp index 5e217c9..e607e33 100644 --- a/examples/flappingBeam_FSI2_coupledFSI_2Dt.cpp +++ b/examples/flappingBeam_FSI2_coupledFSI_2Dt.cpp @@ -84,7 +84,7 @@ int main(int argc, char* argv[]) index_t numUniRef = 3; // time integration real_t timeStep = 0.01; - real_t timeSpan = 15.; + real_t timeSpan = 0.02; //was 15. real_t thetaFluid = 0.5; real_t thetaSolid = 1.; index_t maxCouplingIter = 10; diff --git a/gsBiharmonicAssembler.hpp b/gsBiharmonicAssembler.hpp index 3db2b98..93521ba 100644 --- a/gsBiharmonicAssembler.hpp +++ b/gsBiharmonicAssembler.hpp @@ -17,7 +17,7 @@ #include #include -#include +#include namespace gismo { @@ -89,8 +89,8 @@ void gsBiharmonicAssembler::assemble(bool saveEliminationMatrix) eliminationMatrix.reservePerColumn(m_system.numColNz(m_bases[0],m_options)); } - gsVisitorBiharmonic visitor(*m_pde_ptr, saveEliminationMatrix ? &eliminationMatrix : nullptr); - Base::template push >(visitor); + gsVisitorBiharmonicMixed visitor(*m_pde_ptr, saveEliminationMatrix ? &eliminationMatrix : nullptr); + Base::template push >(visitor); m_system.matrix().makeCompressed(); diff --git a/gsVisitorBiharmonicMixed.h b/gsVisitorBiharmonicMixed.h new file mode 100644 index 0000000..498a360 --- /dev/null +++ b/gsVisitorBiharmonicMixed.h @@ -0,0 +1,170 @@ +/** @file gsVisitorBiharmonicMixed.h + + @brief Visitor class for mixed formulation for the biharmonic equation + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + A.Shamanskiy (2016 - ...., TU Kaiserslautern) +*/ + +#pragma once + +#include +#include + +namespace gismo +{ + +template +class gsVisitorBiharmonicMixed +{ +public: + + gsVisitorBiharmonicMixed(const gsPde & pde_, gsSparseMatrix * elimMatrix = nullptr) + : pde_ptr(static_cast*>(&pde_)), + N_M(), N_A(), localStiffening(), + elimMat(elimMatrix) + {} + + void initialize(const gsBasisRefs & basisRefs, + const index_t patchIndex, + const gsOptionList & options, + gsQuadRule & rule) + { + GISMO_UNUSED(patchIndex); + // a quadrature rule is defined by the basis for the auxiliary variable. + // the same rule is used for the main variable + rule = gsQuadrature::get(basisRefs.back(), options); + // saving necessary info + localStiffening = options.getReal("LocalStiff"); + // resize containers for global indices + globalIndices.resize(2); + blockNumbers.resize(2); + } + + inline void evaluate(const gsBasisRefs & basisRefs, + const gsGeometry & geo, + const gsMatrix & quNodes) + { + // store quadrature points of the element for geometry evaluation + md.points = quNodes; + // NEED_VALUE to get points in the physical domain for evaluation of the RHS + // NEED_MEASURE to get the Jacobian determinant values for integration + // NEED_GRAD_TRANSFORM to get the Jacobian matrix to transform gradient from the parametric to physical domain + md.flags = NEED_VALUE | NEED_MEASURE | NEED_GRAD_TRANSFORM; + // Compute image of the quadrature points plus gradient, jacobian and other necessary data + geo.computeMap(md); + // find local indices of the main and auxiliary basis functions active on the element + basisRefs.front().active_into(quNodes.col(0),localIndicesMain); + N_M = localIndicesMain.rows(); + basisRefs.back().active_into(quNodes.col(0), localIndicesAux); + N_A = localIndicesAux.rows(); + // Evaluate basis functions and their derivatives on the element + basisRefs.front().evalAllDers_into(quNodes,1,basisValuesMain); + basisRefs.back().evalAllDers_into(quNodes,1,basisValuesAux); + // Evaluate right-hand side at the image of the quadrature points + pde_ptr->rhs()->eval_into(md.values[0],forceValues); + } + + inline void assemble(gsDomainIterator & element, + const gsVector & quWeights) + { + GISMO_UNUSED(element); + // Initialize local matrix/rhs // 0 | B^T = L + localMat.setZero(N_M + N_A, N_M + N_A); // --|-- matrix structure + localRhs.setZero(N_M + N_A,pde_ptr->numRhs()); // B | A = 0 + + // Loop over the quadrature nodes + for (index_t q = 0; q < quWeights.rows(); ++q) + { + // Multiply quadrature weight by the geometry measure + const T weightMatrix = quWeights[q] * pow(md.measure(q),1-localStiffening); + const T weightRHS = quWeights[q] * md.measure(q); + // Compute physical gradients of the basis functions at q as a 1 x numActiveFunction matrix + transformGradients(md, q, basisValuesAux[1], physGradAux); + transformGradients(md, q, basisValuesMain[1], physGradMain); + // matrix A + block = weightMatrix * basisValuesAux[0].col(q) * basisValuesAux[0].col(q).transpose(); + localMat.block(N_M,N_M,N_A,N_A) += block.block(0,0,N_A,N_A); + // matrix B + block = weightMatrix * physGradAux.transpose()*physGradMain; // N_A x N_M + localMat.block(0,N_M,N_M,N_A) += block.block(0,0,N_A,N_M).transpose(); + localMat.block(N_M,0,N_A,N_M) += block.block(0,0,N_A,N_M); + // rhs contribution + localRhs.middleRows(0,N_M).noalias() += weightRHS * basisValuesMain[0].col(q) * forceValues.col(q).transpose(); + } + } + + inline void localToGlobal(const int patchIndex, + const std::vector > & eliminatedDofs, + gsSparseSystem & system) + { + // compute global indices + system.mapColIndices(localIndicesMain,patchIndex,globalIndices[0],0); + system.mapColIndices(localIndicesAux,patchIndex,globalIndices[1],1); + blockNumbers.at(0) = 0; + blockNumbers.at(1) = 1; + // push to global system + system.pushToRhs(localRhs,globalIndices,blockNumbers); + system.pushToMatrix(localMat,globalIndices,eliminatedDofs,blockNumbers,blockNumbers); + + // push to the elimination system + if (elimMat != nullptr) + { + index_t globalI,globalElimJ; + index_t elimSize = 0; + for (short_t dJ = 0; dJ < 2; ++dJ) + { + for (short_t dI = 0; dI < 2; ++dI) + for (index_t i = 0; i < N_M; ++i) + if (system.colMapper(dI).is_free_index(globalIndices[dI].at(i))) + { + system.mapToGlobalRowIndex(localIndicesMain.at(i),patchIndex,globalI,dI); + for (index_t j = 0; j < N_M; ++j) + if (!system.colMapper(dJ).is_free_index(globalIndices[dJ].at(j))) + { + globalElimJ = system.colMapper(dJ).global_to_bindex(globalIndices[dJ].at(j)); + elimMat->coeffRef(globalI,elimSize+globalElimJ) += localMat(N_M*dI+i,N_M*dJ+j); + } + } + elimSize += eliminatedDofs[dJ].rows(); + } + } + } + +protected: + // problem info + const gsPoissonPde * pde_ptr; + // geometry mapping + gsMapData md; + // local components of the global linear system + gsMatrix localMat; + gsMatrix localRhs; + // local indices (at the current patch) of basis functions active at the current element + gsMatrix localIndicesMain; + gsMatrix localIndicesAux; + // number of main and auxiliary basis functions active at the current element + index_t N_M, N_A; + // values and derivatives of main basis functions at quadrature points at the current element + // values are stored as a N_M x numQuadPoints matrix; not sure about derivatives, must be smth like N_M x numQuadPoints + // same for the auxiliary basis functions + std::vector > basisValuesMain; + std::vector > basisValuesAux; + // RHS values at quadrature points at the current element; stored as a 1 x numQuadPoints matrix + gsMatrix forceValues; + // all temporary matrices defined here for efficiency + gsMatrix block, physGradMain, physGradAux; + real_t localStiffening; + // elimination matrix to efficiently change Dirichlet degrees of freedom + gsSparseMatrix * elimMat; + // containers for global indices + std::vector< gsMatrix > globalIndices; + gsVector blockNumbers; +}; + +} // namespace gismo diff --git a/gsWriteParaviewMultiPhysics.hpp b/gsWriteParaviewMultiPhysics.hpp index de2dccf..81d534a 100644 --- a/gsWriteParaviewMultiPhysics.hpp +++ b/gsWriteParaviewMultiPhysics.hpp @@ -88,7 +88,7 @@ void gsWriteParaviewMultiPhysics(std::map*> fields fields.begin()->second->igaFunction(i).basis() : fields.begin()->second->patch(i).basis(); gsWriteParaviewMultiPhysicsSinglePatch( fields, i, fn + util::to_string(i), npts); - collection.addPart(baseName + util::to_string(i) + ".vts", -1, "", i ); + collection.addPart(baseName + util::to_string(i) + ".vts", -1, "Solution", i ); if ( mesh ) { @@ -114,7 +114,7 @@ void gsWriteParaviewMultiPhysicsTimeStep(std::map { std::string patchFileName = fn + util::to_string(time) + "_" + util::to_string(p); gsWriteParaviewMultiPhysicsSinglePatch(fields,p,patchFileName,npts); - collection.addPart(gsFileManager::getFilename(patchFileName),time,"solution",p); + collection.addPart(gsFileManager::getFilename(patchFileName),time,"Solution",p); } }