Commit 2c30fc2b authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

removed boost from CMakeLists

parents 281556c2 d53bffff
...@@ -4,10 +4,9 @@ cache: ...@@ -4,10 +4,9 @@ cache:
- install/ - install/
dune:git--gcc: dune:git--gcc:
image: duneci/dune-fufem:git image: mathiwr/dune:latest
script: script:
- ./contrib/ci-setup - dunecontrol --current all
- dunecontrol --opts=/duneci/opts.gcc --current all
- dunecontrol --current make test - dunecontrol --current make test
only: only:
- develop - develop
...@@ -13,6 +13,8 @@ This style-guide is intended for developers writing code for AMDiS, is not compl ...@@ -13,6 +13,8 @@ This style-guide is intended for developers writing code for AMDiS, is not compl
Parts of this convention are taken from well established style guides, like the *Google C++ Style Guide*. Parts of this convention are taken from well established style guides, like the *Google C++ Style Guide*.
In general, the code should follow the [C++ Core Guidelines](https://github.com/isocpp/CppCoreGuidelines).
## File naming conventions ## File naming conventions
Filenames should be mixed lower and upper case, starting with an uppercase letter. They should not include underscores or dashed. Use an uppercase letter to indicate a new subword. Sourcefiles should end in `.cpp` and header files should end in `.hpp`. In case you move the code of a template class to a separate file, included textually at the end of the corresponding header file, use the extensions `.inc.hpp`. Filenames should be mixed lower and upper case, starting with an uppercase letter. They should not include underscores or dashed. Use an uppercase letter to indicate a new subword. Sourcefiles should end in `.cpp` and header files should end in `.hpp`. In case you move the code of a template class to a separate file, included textually at the end of the corresponding header file, use the extensions `.inc.hpp`.
...@@ -26,19 +28,19 @@ Examples of valid filenames: ...@@ -26,19 +28,19 @@ Examples of valid filenames:
* `DOFVector.cpp` * `DOFVector.cpp`
* `DOFVector.inc.hpp` (the implementation of the methods of the template class `DOFVector<T>`) * `DOFVector.inc.hpp` (the implementation of the methods of the template class `DOFVector<T>`)
Do not use filenames that already exist in /usr/include or are stdandard C/C++ include files, such as `math.h` (remember that windows files-systems are case insensitive and thus, there is not difference between `math.h` and `Math.H`.) Do not use filenames that already exist in /usr/include or are stdandard C/C++ include files, such as `math.h` (remember that windows files-systems are case insensitive and thus, there is no difference between `math.h` and `Math.H`.)
## Generale file structure ## Generale file structure
Every header file should start with a copyright notice and an include guard `#pragma once`, where the text of the copyright notice is given in the file `tools/license.templ.txt` and can automatically by added, using the script files in the `tools` directory: Every header file should start with a copyright notice and an include guard `#pragma once`, where the text of the copyright notice is given in the file `tools/license.templ.txt` and can automatically by added, using the script files in the `tools` directory:
``` c++ ``` c++
// Software License for AMDiS2 // Software License for dune-amdis
// //
// Copyright (c) 2015 Institute for Scientific Computing, TU-Dresden // Copyright (c) 2015 Institute for Scientific Computing, Technische Universitaet Dresden
// All rights reserved. // All rights reserved.
// Authors: Simon Praetorius // Authors: Simon Praetorius
// //
// This file is part of the AMDiS2 Library // This file is part of the dune-amdis Library
// see also the LICENSE file in the distribution. // see also the LICENSE file in the distribution.
#pragma once #pragma once
...@@ -55,7 +57,7 @@ All of a project's header files should be listed as descendants of the project's ...@@ -55,7 +57,7 @@ All of a project's header files should be listed as descendants of the project's
* Other external libraries' header files. * Other external libraries' header files.
* Your project's header files. * Your project's header files.
For better readability a comment above each group can be added. Within each section the includes should be ordered alphabetically. Project's header files should be surrounded by `"`, while external header files hould be surrounded by `<...>`. For better readability a comment above each group can be added. Within each section the includes should be ordered alphabetically. Project's header files should be surrounded by `"`, while external header files should be surrounded by `<...>`.
For example, the includes in `io/VtkWriter.cpp` might look like this: For example, the includes in `io/VtkWriter.cpp` might look like this:
...@@ -93,15 +95,15 @@ namespace AMDiS ...@@ -93,15 +95,15 @@ namespace AMDiS
} // end namespace AMDiS } // end namespace AMDiS
``` ```
Implementation details are put into a subnamespace `Impl`. Implementation details are put into a subnamespace `Impl`. A few more subnamespaces of `AMDiS` are allowed, e.g., `Concepts`. If onw of these subnamespaces need another subsubnamespace for implementation details, it should be names `_Impl`.
## Line length ## Line length
Each line of text in your code should be at most 80 characters long. Each line of text in your code should be at most 100 characters long.
**Exceptions:** **Exceptions:**
* An #include statement with a long path may exceed 80 columns. * An #include statement with a long path may exceed 100 columns.
* A raw-string literal may have content that exceeds 80 characters. * A raw-string literal may have content that exceeds 100 characters.
* ... * ...
## Indentation ## Indentation
......
#! /bin/bash
set -e
set -x
root=${PWD}
if [ ! -d install/MTL ]; then
mkdir -p install/MTL
cd install/MTL
curl -o MTL.tar.gz "http://www.simunova.com/downloads/mtl4/MTL-4.0.9555-Linux.tar.gz"
tar --strip-components=2 -xf MTL.tar.gz
rm MTL.tar.gz
cd ${root}
fi
# if [ ! -d install/gtest ]; then
# mkdir -p install/gtest
# cd install/gtest
# cmake -DCMAKE_INSTALL_PREFIX=${root}/install/gtest /usr/src/gtest
# make
# cd ${root}
# fi
...@@ -21,28 +21,27 @@ add_dune_alberta_flags("duneamdis" OBJECT USE_GENERIC) ...@@ -21,28 +21,27 @@ add_dune_alberta_flags("duneamdis" OBJECT USE_GENERIC)
target_compile_definitions("duneamdis" PUBLIC AMDIS_BACKEND_MTL=1) target_compile_definitions("duneamdis" PUBLIC AMDIS_BACKEND_MTL=1)
target_compile_options("duneamdis" PUBLIC -ftemplate-backtrace-limit=0) target_compile_options("duneamdis" PUBLIC -ftemplate-backtrace-limit=0)
set(BOOST_VERSION "1.54") # set(BOOST_VERSION "1.54")
set(BOOST_LIBS_REQUIRED system program_options) # set(BOOST_LIBS_REQUIRED system program_options)
if (NOT BUILD_SHARED_LIBS) # if (NOT BUILD_SHARED_LIBS)
set(Boost_USE_STATIC_LIBS ON) # set(Boost_USE_STATIC_LIBS ON)
endif (NOT BUILD_SHARED_LIBS) # endif (NOT BUILD_SHARED_LIBS)
find_package(Boost ${BOOST_VERSION} REQUIRED ${BOOST_LIBS_REQUIRED}) # find_package(Boost ${BOOST_VERSION} REQUIRED ${BOOST_LIBS_REQUIRED})
if (Boost_FOUND) # if (Boost_FOUND)
add_library(boost INTERFACE) # add_library(boost INTERFACE)
target_include_directories(boost INTERFACE ${Boost_INCLUDE_DIR}) # target_include_directories(boost INTERFACE ${Boost_INCLUDE_DIR})
target_link_libraries(boost INTERFACE ${Boost_LIBRARIES}) # target_link_libraries(boost INTERFACE ${Boost_LIBRARIES})
target_link_libraries(duneamdis INTERFACE boost) # target_link_libraries(duneamdis INTERFACE boost)
if (MSVC_SHARED_LIBS) # if (MSVC_SHARED_LIBS)
link_directories(${Boost_LIBRARY_DIRS}) # link_directories(${Boost_LIBRARY_DIRS})
target_compile_definitions("duneamdis" INTERFACE ${Boost_LIB_DIAGNOSTIC_DEFINITIONS}) # target_compile_definitions("duneamdis" INTERFACE ${Boost_LIB_DIAGNOSTIC_DEFINITIONS})
endif (MSVC_SHARED_LIBS) # endif (MSVC_SHARED_LIBS)
endif (Boost_FOUND) # endif (Boost_FOUND)
find_package(MTL REQUIRED find_package(MTL REQUIRED
PATHS ${CMAKE_SOURCE_DIR}/install/MTL/share/mtl PATHS /usr/local/lib/mtl4)
/iwr/modules/mtl4/share/mtl)
if (MTL_FOUND) if (MTL_FOUND)
target_include_directories("duneamdis" PUBLIC ${MTL_INCLUDE_DIRS}) target_include_directories("duneamdis" PUBLIC ${MTL_INCLUDE_DIRS})
# target_link_libraries("duneamdis" PUBLIC ${MTL_LIBRARIES}) # target_link_libraries("duneamdis" PUBLIC ${MTL_LIBRARIES})
......
...@@ -42,16 +42,16 @@ namespace AMDiS ...@@ -42,16 +42,16 @@ namespace AMDiS
* evaluation of the scalar product `c = C(phi(lambda))*v(lambda)` with lambda a local coordinate. * evaluation of the scalar product `c = C(phi(lambda))*v(lambda)` with lambda a local coordinate.
**/ **/
template <class Term> template <class Term>
Self& addZOT(Term const& c) Operator& addZOT(Term const& c)
{ {
return addZOTImpl(toTerm(c)); return addZOTImpl(toTerm(c));
} }
template <class... Args> template <class... Args>
static shared_ptr<Self> zot(Args&&... args) static Operator zot(Args&&... args)
{ {
auto op = make_shared<Self>(); Self op;
op->addZOT(std::forward<Args>(args)...); op.addZOT(std::forward<Args>(args)...);
return op; return op;
} }
...@@ -61,7 +61,7 @@ namespace AMDiS ...@@ -61,7 +61,7 @@ namespace AMDiS
/// \p type == GRD_PHI and the second one to \p type == GRD_PSI. /// \p type == GRD_PHI and the second one to \p type == GRD_PSI.
/// The coefficient \p b must be a vector expression. /// The coefficient \p b must be a vector expression.
template <class Term> template <class Term>
Self& addFOT(Term const& b, FirstOrderType type) Operator& addFOT(Term const& b, FirstOrderType type)
{ {
return addFOTImpl(toTerm(b), type); return addFOTImpl(toTerm(b), type);
} }
...@@ -71,16 +71,16 @@ namespace AMDiS ...@@ -71,16 +71,16 @@ namespace AMDiS
/// \p type == GRD_PHI and the second one to \p type == GRD_PSI. /// \p type == GRD_PHI and the second one to \p type == GRD_PSI.
/// The coefficient \p b must be a scalar expression. /// The coefficient \p b must be a scalar expression.
template <class Term> template <class Term>
Self& addFOT(Term const& b, std::size_t i, FirstOrderType type) Operator& addFOT(Term const& b, std::size_t i, FirstOrderType type)
{ {
return addFOTImpl(toTerm(b), i, type); return addFOTImpl(toTerm(b), i, type);
} }
template <class... Args> template <class... Args>
static shared_ptr<Self> fot(Args&&... args) static Operator fot(Args&&... args)
{ {
auto op = make_shared<Self>(); Self op;
op->addFOT(std::forward<Args>(args)...); op.addFOT(std::forward<Args>(args)...);
return op; return op;
} }
...@@ -88,7 +88,7 @@ namespace AMDiS ...@@ -88,7 +88,7 @@ namespace AMDiS
/// Add coefficients for second-order operator < grad(psi), A * grad(phi) >, /// Add coefficients for second-order operator < grad(psi), A * grad(phi) >,
/// where \p A can be a matrix or a scalar expression. /// where \p A can be a matrix or a scalar expression.
template <class Term> template <class Term>
Self& addSOT(Term const& A) Operator& addSOT(Term const& A)
{ {
return addSOTImpl(toTerm(A)); return addSOTImpl(toTerm(A));
} }
...@@ -98,16 +98,16 @@ namespace AMDiS ...@@ -98,16 +98,16 @@ namespace AMDiS
/// corresponds to the derivative component of the test-function and the /// corresponds to the derivative component of the test-function and the
/// second index \p _j to the derivative component of the trial function. /// second index \p _j to the derivative component of the trial function.
template <class Term> template <class Term>
Self& addSOT(Term const& A, std::size_t i, std::size_t j) Operator& addSOT(Term const& A, std::size_t i, std::size_t j)
{ {
return addSOTImpl(toTerm(A), i, j); return addSOTImpl(toTerm(A), i, j);
} }
template <class... Args> template <class... Args>
static shared_ptr<Self> sot(Args&&... args) static Operator sot(Args&&... args)
{ {
auto op = make_shared<Self>(); Self op;
op->addSOT(std::forward<Args>(args)...); op.addSOT(std::forward<Args>(args)...);
return op; return op;
} }
......
...@@ -142,10 +142,8 @@ namespace AMDiS ...@@ -142,10 +142,8 @@ namespace AMDiS
void addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j, void addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j,
double* factor = nullptr, double* estFactor = nullptr); double* factor = nullptr, double* estFactor = nullptr);
void addMatrixOperator(std::map< std::pair<int,int>, void addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops);
std::shared_ptr<ElementOperator> > ops); void addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator, bool> > ops);
void addMatrixOperator(std::map< std::pair<int,int>,
std::pair<std::shared_ptr<ElementOperator>, bool> > ops);
/** @} */ /** @} */
...@@ -159,8 +157,8 @@ namespace AMDiS ...@@ -159,8 +157,8 @@ namespace AMDiS
void addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i, void addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i,
double* factor = nullptr, double* estFactor = nullptr); double* factor = nullptr, double* estFactor = nullptr);
void addVectorOperator(std::map< int, std::shared_ptr<ElementOperator> > ops); void addVectorOperator(std::map< int, ElementOperator> ops);
void addVectorOperator(std::map< int, std::pair<std::shared_ptr<ElementOperator>, bool> > ops); void addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops);
/** @} */ /** @} */
/// Adds a Dirichlet boundary condition /// Adds a Dirichlet boundary condition
......
...@@ -121,24 +121,24 @@ addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j, ...@@ -121,24 +121,24 @@ addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j,
template <class Traits> template <class Traits>
void ProblemStat<Traits>:: void ProblemStat<Traits>::
addMatrixOperator(std::map< std::pair<int,int>, std::shared_ptr<ElementOperator> > ops) addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
{ {
for (auto op : ops){ for (auto op : ops){
int r = op.first.first; int r = op.first.first;
int c = op.first.second; int c = op.first.second;
matrixOperators[r][c].push_back({op.second}); matrixOperators[r][c].push_back({std::make_shared<ElementOperator>(op.second)});
matrixChanging[r][c] = true; matrixChanging[r][c] = true;
} }
} }
template <class Traits> template <class Traits>
void ProblemStat<Traits>:: void ProblemStat<Traits>::
addMatrixOperator(std::map< std::pair<int,int>, std::pair<std::shared_ptr<ElementOperator>,bool> > ops) addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops)
{ {
for (auto op : ops) { for (auto op : ops) {
int r = op.first.first; int r = op.first.first;
int c = op.first.second; int c = op.first.second;
matrixOperators[r][c].push_back({op.second.first}); matrixOperators[r][c].push_back({std::make_shared<ElementOperator>(op.second.first)});
matrixChanging[r][c] = matrixChanging[r][c] || op.second.second; matrixChanging[r][c] = matrixChanging[r][c] || op.second.second;
} }
} }
...@@ -163,20 +163,20 @@ addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i, double* ...@@ -163,20 +163,20 @@ addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i, double*
template <class Traits> template <class Traits>
void ProblemStat<Traits>:: void ProblemStat<Traits>::
addVectorOperator(std::map< int, std::shared_ptr<ElementOperator> > ops) addVectorOperator(std::map< int, ElementOperator> ops)
{ {
for (auto op : ops) { for (auto op : ops) {
vectorOperators[op.first].push_back({op.second}); vectorOperators[op.first].push_back({std::make_shared<ElementOperator>(op.second)});
vectorChanging[op.first] = true; vectorChanging[op.first] = true;
} }
} }
template <class Traits> template <class Traits>
void ProblemStat<Traits>:: void ProblemStat<Traits>::
addVectorOperator(std::map< int, std::pair<std::shared_ptr<ElementOperator>, bool> > ops) addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
{ {
for (auto op : ops) { for (auto op : ops) {
vectorOperators[op.first].push_back({op.second.first}); vectorOperators[op.first].push_back({std::make_shared<ElementOperator>(op.second.first)});
vectorChanging[op.first] = vectorChanging[op.first] || op.second.second; vectorChanging[op.first] = vectorChanging[op.first] || op.second.second;
} }
} }
......
...@@ -94,9 +94,9 @@ namespace AMDiS ...@@ -94,9 +94,9 @@ namespace AMDiS
solverMpL2->solve(MpL2, x, b, solverInfo); solverMpL2->solve(MpL2, x, b, solverInfo);
} }
MTLMatrix const& getM() const { return matM ? *matM : (*Super::A)(2_c, 2_c); } // < u, v > MTLMatrix const& getM() const { return matM ? *matM : (*Super::A)(_2, _2); } // < u, v >
MTLMatrix const& getL0() const { return matL0 ? *matL0 : (*Super::A)(1_c, 0_c); } // < M0*tau*grad(u), grad(v) > MTLMatrix const& getL0() const { return matL0 ? *matL0 : (*Super::A)(_1, _0); } // < M0*tau*grad(u), grad(v) >
MTLMatrix const& getL() const { return matL ? *matL : (*Super::A)(2_c, 1_c); } // < grad(u), grad(v) > MTLMatrix const& getL() const { return matL ? *matL : (*Super::A)(_2, _1); } // < grad(u), grad(v) >
double getTau() const { return *tau; } double getTau() const { return *tau; }
......
...@@ -13,9 +13,9 @@ ...@@ -13,9 +13,9 @@
using namespace AMDiS; using namespace AMDiS;
// 1 component with polynomial degree 1 // 1 component with polynomial degree 1
using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
using ElliptParam = LagrangeTraits<Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>, 1>; using ElliptParam = LagrangeTraits<Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>, 1>;
using ElliptProblem = ProblemStat<ElliptParam>; using ElliptProblem = ProblemStat<ElliptParam>;
using ElliptProblemInstat = ProblemInstat<ElliptParam>;
int main(int argc, char** argv) int main(int argc, char** argv)
{ {
...@@ -24,9 +24,6 @@ int main(int argc, char** argv) ...@@ -24,9 +24,6 @@ int main(int argc, char** argv)
ElliptProblem prob("ellipt"); ElliptProblem prob("ellipt");
prob.initialize(INIT_ALL); prob.initialize(INIT_ALL);
ElliptProblemInstat probInstat("ellipt", prob);
probInstat.initialize(INIT_UH_OLD);
AdaptInfo adaptInfo("adapt"); AdaptInfo adaptInfo("adapt");
using Op = ElliptProblem::ElementOperator; using Op = ElliptProblem::ElementOperator;
...@@ -35,7 +32,7 @@ int main(int argc, char** argv) ...@@ -35,7 +32,7 @@ int main(int argc, char** argv)
opL.addSOT(1.0); opL.addSOT(1.0);
prob.addMatrixOperator(opL, 0, 0); prob.addMatrixOperator(opL, 0, 0);
opForce.addZOT( eval([](auto const& x) { return -1.0; }) ); opForce.addZOT([](auto const& x) { return -1.0; });
prob.addVectorOperator(opForce, 0); prob.addVectorOperator(opForce, 0);
......
...@@ -25,13 +25,12 @@ int main(int argc, char** argv) ...@@ -25,13 +25,12 @@ int main(int argc, char** argv)
{ {
AMDiS::init(argc, argv); AMDiS::init(argc, argv);
// Create grid from structured grid // Create grid from structured grid
std::array<unsigned int, 2> n = {{4, 4}}; std::array<unsigned int, 2> n = {{4, 4}};
Dune::FieldVector<double, 2> lower = {{0.0, 0.0}}; Dune::FieldVector<double, 2> lower = {{0.0, 0.0}};
Dune::FieldVector<double, 2> upper = {{1.0, 1.0}}; Dune::FieldVector<double, 2> upper = {{1.0, 1.0}};
auto grid = Dune::StructuredGridFactory<Grid>::createSimplexGrid(lower, upper, n); auto grid = Dune::StructuredGridFactory<Grid>::createSimplexGrid(lower, upper, n);
// NOTE: can not be used with AlbertaGrid
ElliptProblem prob("ellipt", *grid); ElliptProblem prob("ellipt", *grid);
prob.initialize(INIT_ALL); prob.initialize(INIT_ALL);
...@@ -63,13 +62,6 @@ int main(int argc, char** argv) ...@@ -63,13 +62,6 @@ int main(int argc, char** argv)
*prob.getSolution() = 0.0; // maybe not necessary *prob.getSolution() = 0.0; // maybe not necessary
prob.buildAfterCoarsen(adaptInfo, Flag(0)); prob.buildAfterCoarsen(adaptInfo, Flag(0));
// write matrix to file
{
mtl::io::matrix_market_ostream out("matrix.mtx");
out << prob.getSystemMatrix()->getMatrix<0,0>();
}
prob.solve(adaptInfo); prob.solve(adaptInfo);
prob.writeFiles(adaptInfo, true); prob.writeFiles(adaptInfo, true);
......
...@@ -58,10 +58,10 @@ int main(int argc, char** argv) ...@@ -58,10 +58,10 @@ int main(int argc, char** argv)
opL.addSOT( constant(1.0) ); opL.addSOT( constant(1.0) );
prob.addMatrixOperator(opL, 0, 0); prob.addMatrixOperator(opL, 0, 0);
opTimeRhs.addZOT( valueOf(prob.getSolution(0_c)) ); opTimeRhs.addZOT( valueOf(prob.getSolution(_0)) );
prob.addVectorOperator(opTimeRhs, 0, probInstat.getInvTau()); prob.addVectorOperator(opTimeRhs, 0, probInstat.getInvTau());
opForce.addZOT( eval([](auto const& x) { return -1.0; }) ); opForce.addZOT([](auto const& x) { return -1.0; });
prob.addVectorOperator(opForce, 0); prob.addVectorOperator(opForce, 0);
......
...@@ -50,38 +50,39 @@ int main(int argc, char** argv) ...@@ -50,38 +50,39 @@ int main(int argc, char** argv)
using Op = StokesProblem::ElementOperator; using Op = StokesProblem::ElementOperator;
forEach(range_<0,DOW>, [&](auto const _i) forEach(range_<0,DOW>, [&](auto const _i)
{ {
// <1/tau * u_i, v_i> // <1/tau * u_i, v_i>
auto opTime = Op::zot( density ); auto opTime = Op::zot( density );
auto opTimeOld = Op::zot( valueOf(prob.getSolution(_i), density) ); auto opTimeOld = Op::zot( valueOf(prob.getSolution(_i), density) );
prob.addMatrixOperator(*opTime, _i, _i, probInstat.getInvTau()); prob.addMatrixOperator(opTime, _i, _i, probInstat.getInvTau());
prob.addVectorOperator(*opTimeOld, _i, probInstat.getInvTau()); prob.addVectorOperator(opTimeOld, _i, probInstat.getInvTau());
#if STRATEGY == 1
// <(u * nabla)u_i^old, v_i> // <(u * nabla)u_i^old, v_i>
forEach(range_<0, DOW>, [&](auto const _j) forEach(range_<0, DOW>, [&](auto const _j)
{ {
auto opNonlin = Op::zot( derivativeOf(prob.getSolution(_i), _j, density) ); auto opNonlin = Op::zot( derivativeOf(prob.getSolution(_i), _j, density) );
prob.addMatrixOperator(*opNonlin, _i, _j); prob.addMatrixOperator(opNonlin, _i, _j);
}); });
#elif STRATEGY == 2
// <(u^old * nabla)u_i, v_i> // <(u^old * nabla)u_i, v_i>
forEach(range_<0, DOW>, [&](auto const _j) forEach(range_<0, DOW>, [&](auto const _j)
{ {
auto opNonlin = Op::fot( valueOf(prob.getSolution(_j), density), _j, GRD_PHI ); auto opNonlin = Op::fot( valueOf(prob.getSolution(_j), density), _j, GRD_PHI );
prob.addMatrixOperator(*opNonlin, _i, _i); prob.addMatrixOperator(opNonlin, _i, _i);
}); });
#else
// <(u^old * nabla)u_i^old, v_i> // <(u^old * nabla)u_i^old, v_i>
forEach(range_<0, DOW>, [&](auto const _j) forEach(range_<0, DOW>, [&](auto const _j)
{ {
auto opNonlin = Op::zot( density * valueOf(prob.getSolution(_j)) * derivativeOf(prob.getSolution(_i), _j) ); auto opNonlin = Op::zot( func([density](double uj, double dj_ui)
prob.addVectorOperator(*opNonlin, _i); {
return density * uj * dj_ui;
}, valueOf(prob.getSolution(_j)), derivativeOf(prob.getSolution(_i), _j)) );
prob.addVectorOperator(opNonlin, _i);
}); });
#endif
auto opL = std::make_pair( Op::sot( viscosity ), false ); // <viscosity*grad(u_i), grad(v_i)> auto opL = std::make_pair( Op::sot( viscosity ), false ); // <viscosity*grad(u_i), grad(v_i)>
auto opP = std::make_pair( Op::fot( 1.0, _i, GRD_PSI ), false); // <p, d_i(v_i)>