Commit fa33375e authored by Lukas Riedel's avatar Lukas Riedel 🎧
Browse files

observation operator now works for single positions

parent 3a7753be
......@@ -166,44 +166,46 @@ public:
return res;
}
void observation_operator ()
std::vector<RF> observation_operator (const typename Traits::Domain pos)
{
SolutionContainer kfu_raw = Dune::PDELab::Backend::native(kfu);
SolutionContainer kfu_raw = cache_water_content();
std::vector<RF> res(kfu_raw.size());
typedef LocalFunctionSpace<KFGFS> LFS;
typedef LFSIndexCache<LFS> LFSCache;
typedef typename KFU::template ConstLocalView<LFSCache> XView;
typedef FiniteElementInterfaceSwitch<
typename Dune::PDELab::LocalFunctionSpace<KFGFS>::Traits::FiniteElementType
> FESwitch;
LFS lfs(kfgfs);
LFSCache lfs_cache(lfs);
XView x_view(kfu);
std::vector<RF> xl(kfgfs.maxLocalSize());
//std::vector<typename Traits::RangeType> yb(this->gfs->maxLocalSize());
typedef FiniteElementInterfaceSwitch<
typename Dune::PDELab::LocalFunctionSpace<KFGFS>::Traits::FiniteElementType
> FESwitch;
for(auto it = kfgv.template begin<0>(); it != kfgv.template end<0>(); ++it){
auto e = *it;
lfs.bind(e);
lfs_cache.update();
x_view.bind(lfs_cache);
x_view.read(xl);
lfs.debug();
for(std::size_t i = 0; i<xl.size(); ++i){
//if(xl[i] == kfu[lfs_cache.containerIndex(i)]){
if(xl[i] == x_view[i]){
std::cout << "local vector index " << i << " maps to global index " << x_view.cache().containerIndex(i)[0] << std::endl;
}
else{
std::cout << "-> local vector index " << i << " and global index " << x_view.cache().containerIndex(i)[0] << " don't match!" << std::endl;
}
auto global_index = x_view.cache().containerIndex(i)[0];
std::cout << "xl["<<i<<"] = " << xl[i] << ", " << "kfu_raw["<<global_index<<"] = " << kfu_raw[global_index] << std::endl;
}
std::vector<typename Traits::Scalar> fem_factors(kfgfs.maxLocalSize()); // FEM coefficient container
// find entity containing the coordinate
auto it = kfgv.template begin<0>();
auto pos_local = pos;
for( ; it != kfgv.template end<0>(); ++it){
const auto geo_type = it->type();
const auto& ref = Dune::ReferenceElements<typename Traits::DF,Traits::dim>::general(geo_type);
pos_local = it->geometry().local(pos);
if(ref.checkInside(pos_local))
break;
}
if(it == kfgv.template end<0>())
DUNE_THROW(Dune::Exception,"Measurement position is not inside grid domain!");
lfs.bind(*it);
lfs_cache.update();
x_view.bind(lfs_cache);
FESwitch::basis(lfs.finiteElement()).evaluateFunction(pos_local,fem_factors);
for(std::size_t i = 0; i<fem_factors.size(); ++i){
const auto index = x_view.cache().containerIndex(i)[0];
res[index] = fem_factors[i][0];
}
return res;
}
};
......
......@@ -36,15 +36,33 @@ int main(int argc, char **argv)
kfs.solution_to_water_content();
auto start = kfs.cache_water_content();
kfs.observation_operator();
Dune::FieldVector<double,2> pos1({0.5,1.8});
Dune::FieldVector<double,2> pos2({0.5,0.2});
auto obs = kfs.observation_operator(pos2);
obs = kfs.observation_operator(pos1);
double val = .0;
for(std::size_t i=0; i<obs.size(); ++i){
val += start[i]*obs[i];
}
std::cout << "value = " << val << std::endl;
kfs.write_vtk_output();
for(int i=0; i<5; i++){
for(int i=0; i<20; i++){
kfs.compute_time_step();
kfs.write_vtk_output();
kfs.solution_to_water_content();
start = kfs.cache_water_content();
double val = .0;
for(std::size_t i=0; i<obs.size(); ++i){
val += start[i]*obs[i];
}
std::cout << "value = " << val << std::endl;
}
auto stop = kfs.propagate(start,1.0E4,2.0E4);
//auto stop = kfs.propagate(start,1.0E4,2.0E4);
//Dune::Dorie::SimpleSimulation<Dune::Dorie::BaseTraits<YaspGrid,Dune::GeometryType::BasicType::cube,2,1>> sim1(helper,yasp_grid,inifile);
//sim1.run();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment