diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8aa39abbedd6b7a15357b61d7b0dc586a4d56e0e..7bc72217e2143b8cf5e795459c3cd5b2b30657b7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -185,6 +185,15 @@ endif(NOT GMSH_INC)
list(APPEND EXTRA_INCS ${GMSH_INC})
list(APPEND EXTRA_LIBS ${GMSH_LIB})
+# Yaml-cpp
+option(ENABLE_YAML_CPP "Enable YAML-CPP library" ON)
+
+if(ENABLE_YAML_CPP)
+ find_package(yaml-cpp REQUIRED)
+ add_definitions(-DENABLE_YAML_CPP)
+endif()
+
+
# GmshFem
find_library(GMSHFEM_LIB gmshfem)
@@ -211,10 +220,15 @@ find_package(OpenMP)
if(OPENMP_FOUND)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
endif(OPENMP_FOUND)
-
if(APPLE AND EXISTS "/usr/local/opt/libomp")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Xpreprocessor -fopenmp -I/usr/local/opt/libomp/include")
list(APPEND EXTRA_LIBS "-L/usr/local/opt/libomp/lib -lomp")
+elseif(APPLE AND EXISTS "/opt/local/lib/libomp")
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Xpreprocessor -fopenmp -I/opt/local/include/libomp")
+ list(APPEND EXTRA_LIBS "-L/opt/local/lib/libomp -lomp")
+elseif(APPLE AND EXISTS "/opt/homebrew/opt/libomp")
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Xpreprocessor -fopenmp -I/opt/homebrew/opt/libomp/include")
+ list(APPEND EXTRA_LIBS "-L/opt/homebrew/opt/libomp/lib -lomp")
endif()
include_directories(${EXTRA_INCS})
@@ -222,28 +236,38 @@ include_directories(${EXTRA_INCS})
file( GLOB LIB_SOURCES contrib/blossom5v2.05/*.cpp contrib/blossom5v2.05/*/*.cpp common/*/*/*.cpp common/*/*.cpp common/*.cpp specific/*/*/*.cpp specific/*/*.cpp specific/*.cpp)
file( GLOB LIB_HEADERS contrib/blossom5v2.05/*.h contrib/blossom5v2.05/*/*.h common/*/*/*.h common/*/*.h common/*.h specific/*/*/*.h specific/*/*.h specific/*.h)
add_library( fwi ${LIB_SOURCES} ${LIB_HEADERS})
+target_link_libraries(fwi yaml-cpp)
-add_executable( synthetics synthetics.cpp ${EXTRA_INCS})
+add_executable( synthetics synthetics.cpp )
target_link_libraries( synthetics fwi ${EXTRA_LIBS})
-add_executable( time_synthetics_wave time_synthetics_wave.cpp ${EXTRA_INCS})
-target_link_libraries( time_synthetics_wave fwi ${EXTRA_LIBS})
-add_executable( directional directional.cpp ${EXTRA_INCS})
+#add_executable( time_synthetics_wave time_synthetics_wave.cpp ${EXTRA_INCS})
+#target_link_libraries( time_synthetics_wave fwi ${EXTRA_LIBS})
+add_executable( directional directional.cpp )
target_link_libraries( directional fwi ${EXTRA_LIBS})
-add_executable( directional_field directional_field.cpp ${EXTRA_INCS})
+add_executable( directional_field directional_field.cpp )
target_link_libraries( directional_field fwi ${EXTRA_LIBS})
-add_executable( gradient gradient.cpp ${EXTRA_INCS})
+add_executable( gradient gradient.cpp )
target_link_libraries( gradient fwi ${EXTRA_LIBS})
#add_executable( convergence convergence.cpp ${EXTRA_INCS})
#target_link_libraries( convergence fwi ${EXTRA_LIBS})
#add_executable( preconditioner preconditioner.cpp ${EXTRA_INCS})
#target_link_libraries( preconditioner fwi ${EXTRA_LIBS})
-add_executable( ip_comparison ip_comparison.cpp ${EXTRA_INCS})
+add_executable( ip_comparison ip_comparison.cpp )
target_link_libraries( ip_comparison fwi ${EXTRA_LIBS})
#add_executable( ob_comparison ob_comparison.cpp ${EXTRA_INCS})
#target_link_libraries( ob_comparison fwi ${EXTRA_LIBS})
-add_executable( multiscale multiscale.cpp ${EXTRA_INCS})
+add_executable( multiscale multiscale.cpp )
target_link_libraries( multiscale fwi ${EXTRA_LIBS})
-add_executable( inversion inversion.cpp ${EXTRA_INCS})
+add_executable( inversion inversion.cpp )
target_link_libraries( inversion fwi ${EXTRA_LIBS})
-add_executable( error error.cpp ${EXTRA_INCS})
+add_executable( error error.cpp)
target_link_libraries( error fwi ${EXTRA_LIBS})
+
+
+# Unit tests
+
+find_package(Boost 1.67.0 COMPONENTS unit_test_framework REQUIRED)
+include_directories(${Boost_INCLUDE_DIRS})
+
+add_executable(tests tests/test.cpp ${EXTRA_INCS})
+target_link_libraries(tests ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} ${EXTRA_LIBS} fwi)
diff --git a/common/configuration.cpp b/common/configuration.cpp
index 949e9dc77de9542177b68d4e2271aa212091b7c0..c8b5ba2b06679ba15849e1ad36f8cda0618b8478 100644
--- a/common/configuration.cpp
+++ b/common/configuration.cpp
@@ -6,11 +6,36 @@
//GmshFWI Library
#include "configuration.h"
+#include <gmsh.h>
+
+namespace gmodel = gmsh::model;
+namespace factory = gmsh::model::geo;
using namespace gmshfem;
using namespace gmshfem::common;
using namespace gmshfem::domain;
+void ConfigurationInterface::mesh() const
+{
+ msg::print << "Generate meshes" << msg::endl;
+ gmsh::option::setNumber("General.Terminal", 1);
+
+ gmodel::add(_name);
+ if (!_remesh)
+ {
+ gmsh::open(_name + ".msh");
+ return;
+ }
+
+ wave_mesh();
+ data_mesh();
+
+ factory::synchronize();
+ gmodel::mesh::generate();
+
+ gmsh::write(_name + ".msh");
+}
+
bool ConfigurationInterface::recIsValid(unsigned int shot,unsigned int rec) const
{
if( !(shot<_ns) && !(rec<nr(shot)) )
@@ -31,7 +56,7 @@ bool ConfigurationInterface::pntIsValid(unsigned int e_r) const
{
if( !(e_r<_np) )
{
- throw Exception("Emitter-receiver "+std::to_string(e_r)+" is out of scope.");
+ throw Exception("Emitter-receiver "+std::to_string(e_r)+" is out of scope. There are only "+std::to_string(_np)+" points.");
}
return true;
}
diff --git a/common/configuration.h b/common/configuration.h
index 1c883eef0ca90f7d5cf5da2ae346a7042758d1cd..0bc01eb06fbb3513a2fe08b1098a184eb4ebb6b3 100644
--- a/common/configuration.h
+++ b/common/configuration.h
@@ -9,6 +9,7 @@
#include "Domain.h"
//Standard Library
#include <array>
+#include <map>
//GmshFWI Library
#include "model/element.h"
#include "model/parametrization.h"
@@ -24,6 +25,7 @@ protected:
std::array<gmshfem::domain::Domain,2> _wave_omega;
std::array<gmshfem::domain::Domain,2> _model_known;
std::array<gmshfem::domain::Domain,2> _model_unknown;
+ std::map<std::string, gmshfem::domain::Domain> _named_domains;
gmshfem::domain::Domain _data_omega;
unsigned int _np;
@@ -52,9 +54,15 @@ public:
gmshFem.userDefinedParameter(remesh, "remesh");
_remesh = ((bool) remesh);
};
+
+ ConfigurationInterface() : _name("dummy"), _parametrization(nullptr), _mc(1,0.)
+ {
+ _remesh = true;
+ };
+
virtual ~ConfigurationInterface() {delete _parametrization;};
- virtual void mesh() const = 0;
+ virtual void mesh() const;
virtual double area() const = 0;
virtual double data_area() const {return gmshfem::post::integrate(1._d_sf,_data_omega,"Gauss0");};
virtual double datapoint_area() const = 0;
@@ -77,6 +85,7 @@ public:
gmshfem::domain::Domain model_evaldomain() const {return _model_unknown[Support::BLK] | _model_unknown[Support::BND];};
gmshfem::domain::Domain wave_evaldomain() const {return _wave_omega[Support::BLK] | _wave_omega[Support::BND] | _points;};
gmshfem::domain::Domain data_evaldomain() const {return _data_omega;};
+ const std::map<std::string, gmshfem::domain::Domain>& named_domains() const {return _named_domains;}
gmshfem::domain::Domain emitter(unsigned int shot, unsigned int em) const;
gmshfem::domain::Domain receiver(unsigned int shot, unsigned int rec) const;
diff --git a/common/data/element.h b/common/data/element.h
index 68be381fe302d8e67c5b6e50cd500622d6ca7071..d17c7c3ad16809526f44c63048a7abcda8073bfa 100644
--- a/common/data/element.h
+++ b/common/data/element.h
@@ -118,6 +118,45 @@ std::complex<double> l2scalp(const Data<T_Physic>& d1,const Data<T_Physic>& d2)
}
return scalp;
}
+template<Physic T_Physic>
+std::complex<double> l2scalpNoConj(const Data<T_Physic>& d1,const Data<T_Physic>& d2)
+{
+ areCompatible(d1,d2);
+ std::complex<double> scalp = 0.;
+ for (unsigned int f = 0; f < d1.nf(); f++)
+ {
+ for (unsigned int s = 0; s < d1.ns(); s++)
+ {
+ for (unsigned int r = 0; r < d1.nr(s); r++)
+ {
+ scalp += d1.value(f,s,r) * d2.value(f,s,r);
+ }
+ }
+ }
+ return scalp;
+}
+// Inner product are not feasible for that physic
+template<>
+inline std::complex<double> l2scalp<Physic::electromagnetic>(const Data<Physic::electromagnetic>& d1,const Data<Physic::electromagnetic>& d2)
+{
+ throw gmshfem::common::Exception("Electromagnetism isn't supported here");
+}
+template<>
+inline std::complex<double> l2scalp<Physic::elastodynamic>(const Data<Physic::elastodynamic>& d1, const Data<Physic::elastodynamic>& d2)
+{
+ throw gmshfem::common::Exception("Elastodynamics isn't supported here");
+}
+template<>
+inline std::complex<double> l2scalpNoConj<Physic::electromagnetic>(const Data<Physic::electromagnetic>& d1,const Data<Physic::electromagnetic>& d2)
+{
+ throw gmshfem::common::Exception("Electromagnetism isn't supported here");
+}
+template<>
+inline std::complex<double> l2scalpNoConj<Physic::elastodynamic>(const Data<Physic::elastodynamic>& d1, const Data<Physic::elastodynamic>& d2)
+{
+ throw gmshfem::common::Exception("Elastodynamics isn't supported here");
+}
+
template<Physic T_Physic>
Data<T_Physic> conj(const Data<T_Physic>& d)
{
diff --git a/common/model/parametrization.h b/common/model/parametrization.h
index 6813da8fb9d521b1aa5339b9c10fd3d7b15dbffd..ed84ae68d93af22fe3970aafaa051f0bd1a1423b 100644
--- a/common/model/parametrization.h
+++ b/common/model/parametrization.h
@@ -52,14 +52,14 @@ public:
return out;
};
- /* operator() (Only perturbed forward): from external to internal, ms in external parametrization */
+ /* operator() (Only perturbed forward): from external to internal, ms in external parametrization. fct is the perturbation */
virtual ModelFunction operator()(const ModelFunction& fct, const ModelFunction& m) const = 0;
ModelFunction operator()(const ModelField& field, const ModelFunction& m) const
{
return operator()(ModelFunction(field),m);
};
- /* operator() (Only adjoint and perturbed adjoint): from internal to external, ms is internal parametrization */
+ /* operator() (Only adjoint and perturbed adjoint): from internal to external, ms is internal parametrization. Fct contains derivatives in internal. */
virtual Sensitivity operator()(const SensitivityStateEvaluator& fct, Order order, Support support, const ModelStateEvaluator& ms) const = 0;
/* size */
diff --git a/common/model/updater.cpp b/common/model/updater.cpp
index eaff8e63acde1cea9c0b984e5af2e44d496fb072..3c904bb70b743f29da5068be93c020a9ab607535 100644
--- a/common/model/updater.cpp
+++ b/common/model/updater.cpp
@@ -70,6 +70,8 @@ void ModelUpdater::update(Type type, const ModelFunction* const m, const ModelFu
void ModelUpdater::update(Type type, const ModelField* const m , const ModelField* const dm )
{
if(_ms._isUpToDate[type]){return;}
+ assert(_innerproduct);
+ assert(_su);
switch (type)
{
@@ -84,11 +86,13 @@ void ModelUpdater::update(Type type, const ModelField* const m , const ModelFiel
}
break;
case Type::AS:
- update(Type::FS,m,dm);
- _ms._field[type] = _innerproduct->update(
- _su->get(to_order(type),Support::BLK,_ms),
- _su->get(to_order(type),Support::BND,_ms)
- );
+ {
+ update(Type::FS, m, dm);
+ auto &bulkSensitivity = _su->get(to_order(type), Support::BLK, _ms);
+ auto &boundarySensitivity = _su->get(to_order(type), Support::BND, _ms);
+
+ _ms._field[type] = _innerproduct->update(bulkSensitivity, boundarySensitivity);
+ }
break;
case Type::PFS:
if(dm==nullptr)
diff --git a/common/statefunctional.cpp b/common/statefunctional.cpp
index 5ce04598d83f8cdc3a68607ac0ae09ec72797307..b519bfa6640a12dd6a3da48b1ead6d6be4061d72 100644
--- a/common/statefunctional.cpp
+++ b/common/statefunctional.cpp
@@ -227,6 +227,7 @@ template<Physic T_Physic>
void StateFunctional<T_Physic>::setModelPerturbation(const ModelField& dm)
{
modelPerturbationIsObsolete();
+ // Update model PFS, e.g. copy the field to the model perturbation object
std::array<bool,4> ToBeUpdated = {false,true,false,false};
_mu.get(ToBeUpdated,nullptr,&dm);
}
@@ -304,6 +305,8 @@ double StateFunctional<T_Physic>::performance_regularization()
template<Physic T_Physic>
const ModelField& StateFunctional<T_Physic>::gradient()
{
+ // Update the model adjoint state and return it.
+ // TODO: no model sent but FS is updated -> Is it always up to date ?
std::array<bool,4> ToBeUpdated = {false,false,true,false};
return _mu.get(ToBeUpdated).field(Type::AS);
}
@@ -378,6 +381,21 @@ double StateFunctional<T_Physic>::directional1(const ModelField &dm2)
return std::real( _innerproduct->product(gradient(),dm2) );
}
+template <Physic T_Physic>
+double StateFunctional<T_Physic>::directional1direct(const ModelField &dm)
+{
+ setModelPerturbation(dm);
+ // Needs data AS (gradient kernel of the objective) and data PFS (perturbation on the receivers)
+ std::array<bool, 5> DataToBeUpdated = {false, true, true, false, false};
+ std::array<bool, 4> ModelToBeUpdated = {false, false, false, false};
+
+ auto &du = _du.get(DataToBeUpdated, _mu.get(ModelToBeUpdated));
+ auto &dataAS = du.state(Type::AS);
+ auto &dataPFS = du.state(Type::PFS);
+
+ return real(l2scalpNoConj(dataAS, dataPFS));
+}
+
/*
* Second order
*/
diff --git a/common/statefunctional.h b/common/statefunctional.h
index 267cfdd858333727b22d30537147595489ae3bea..bc3aea3169afbdffabe67ef7848c5bf721d5bfbb 100644
--- a/common/statefunctional.h
+++ b/common/statefunctional.h
@@ -79,6 +79,8 @@ public:
/* directional1 */
virtual double directional1(const ModelFunction &dm);
virtual double directional1(const ModelField &dm);
+ // Solve the PFS for a direction and compute the variation of cost.
+ virtual double directional1direct(const ModelField &dm);
/*
* Second order
diff --git a/common/wave/element.cpp b/common/wave/element.cpp
index 15a818561b51d37f2f5a62a2709bd5d968a021a1..d074e2c14bcdfedb990376662171b5661c8d9e3b 100644
--- a/common/wave/element.cpp
+++ b/common/wave/element.cpp
@@ -11,9 +11,9 @@ template<Physic T_Physic>
void WaveMultiField<T_Physic>::initializeField(const WaveField<T_Physic>& field)
{
_field = field;
- for (auto it = _dofValues.begin(); it != _dofValues.end(); it++)
+ for (unsigned i = 0; i < _dofValues.size(); ++i)
{
- it->resize(_field.numberOfDofs());
+ _dofValues[i].resize(_field.numberOfDofs());
}
_isInitialized = true;
}
@@ -52,19 +52,29 @@ template<Physic T_Physic>
void WaveAuxilaryField<T_Physic>::assignValues(const std::vector< std::complex<double> > &values)
{
if(_systemValuesToDofValues.size() == 0) {
- _systemValuesToDofValues.reserve(values.size());
+ _systemValuesToDofValues.reserve(_wmf->_field.numberOfDofs());
+ _unknownIndices.reserve(values.size());
+ unsigned counter = 0;
for(auto it = _wmf->_field.begin(); it != _wmf->_field.end(); it++)
{
- if(it->first->type() == gmshfem::dofs::Type::Unknown)
- {
- _systemValuesToDofValues.push_back(it->first->numDof() - 1);
- }
+ // Save all DOFs, including the fixed ones, but remember to only update the unknown ones
+ _systemValuesToDofValues.push_back(it->first->numDof() - 1);
+ if (it->first->type() == gmshfem::dofs::Type::Unknown) {
+ _unknownIndices.push_back(counter);
+ }
+ ++counter;
+
}
}
- for(unsigned int i = 0; i < values.size(); ++i) {
+
+ for (unsigned int i : _unknownIndices)
+ {
_wmf->_dofValues[_index][i] = values[_systemValuesToDofValues[i]];
}
- if(_index == _wmf->_index){_wmf->_isUpToDate=false;}
+ if (_index == _wmf->_index)
+ {
+ _wmf->_isUpToDate = false;
+ }
}
template class WaveMultiField<Physic::acoustic>;
diff --git a/common/wave/element.h b/common/wave/element.h
index 26731b9191f34058061789c66dc12ea9b13d1862..23b38cec2b805847be35981099dfd1a34621ab10 100644
--- a/common/wave/element.h
+++ b/common/wave/element.h
@@ -52,6 +52,7 @@ private:
WaveMultiField<T_Physic>* _wmf;
unsigned int _index;
std::vector< unsigned int > _systemValuesToDofValues;
+ std::vector< unsigned int > _unknownIndices;
public:
WaveAuxilaryField(std::string name, gmshfem::domain::Domain domain, std::string gmodel, const wave::Discretization<T_Physic>& w_discret): WaveField<T_Physic>(name,domain,gmodel,w_discret), _wmf(nullptr), _index(0), _systemValuesToDofValues() {};
void setField(WaveMultiField<T_Physic>* wmf){_wmf=wmf;_systemValuesToDofValues.clear();};
diff --git a/common/wave/equation/equation.cpp b/common/wave/equation/equation.cpp
index 780b1a30b33d0770cfb9530db82d0b1efe3abbed..6f1934a60af0015f35a5d4f0acdacd79e27f5edf 100644
--- a/common/wave/equation/equation.cpp
+++ b/common/wave/equation/equation.cpp
@@ -40,13 +40,13 @@ EquationInterface<T_Physic>::EquationInterface(const ConfigurationInterface* con
{
throw Exception("Equation integration degree (boundary) could not be found.");
}
- double buffer=0.;
+ int buffer=0;
if(!gmshFem.userDefinedParameter(buffer, "equation_boundary"))
{
msg::warning << "Boundary are not considered in equations (default behaviour)." << msg::endl;
}
- if(buffer==0.){_boundary=false;}
- else{_boundary=true;}
+ if(buffer==0){_boundary=false; msg::warning << "Boundary are not enabled in equations." << msg::endl;}
+ else{_boundary=true; msg::warning << "Boundary are enabled in equations." << msg::endl;}
};
template<Physic T_Physic>
diff --git a/convergence.cpp b/convergence.cpp
index 7b08d37a9584797a1a544e49302a7519db15d3d2..7b68f25843dca8f471785aae0493fd3ce1741482 100644
--- a/convergence.cpp
+++ b/convergence.cpp
@@ -41,7 +41,7 @@ int convergence(const GmshFem& gmshFem)
ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
- ConfigurationInterface* configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
std::vector<double> frequency(1,0.);
if(!gmshFem.userDefinedParameter(frequency[0], "frequency0"))
@@ -50,15 +50,15 @@ int convergence(const GmshFem& gmshFem)
}
std::vector<std::string> filename (1,"noname");
gmshFem.userDefinedParameter(filename[0], "data0");
- const Data<T_Physic> d0(filename,frequency,configuration);
+ const Data<T_Physic> d0(filename,frequency,configuration.get());
- ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(d0,configuration,gmshFem);
+ ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(d0,configuration.get(),gmshFem);
wave::Discretization<T_Physic> w_discret(gmshFem);
- std::vector<EquationInterface<T_Physic>*> pequation(1,new ParametrizedEquation<T_Physic>(parametrization,0,2.*M_PI*frequency[0],configuration,w_discret,gmshFem));
+ std::vector<EquationInterface<T_Physic>*> pequation(1,new ParametrizedEquation<T_Physic>(parametrization,0,2.*M_PI*frequency[0],configuration.get(),w_discret,gmshFem));
- StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,nullptr,frequency,pequation,objective);
+ StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration.get(),nullptr,frequency,pequation,objective);
FunctionalInterface* const functional = statefunctional;
double eps0,epsN;
@@ -156,7 +156,6 @@ int convergence(const GmshFem& gmshFem)
delete statefunctional;
delete objective;
delete pequation[0];
- delete configuration;
delete parametrization;
return 0;
diff --git a/directional.cpp b/directional.cpp
index a8a6cd26bc35605d8e4e9ec8e1273ae48848b3e7..53e3771ec05dcb5f94a13e77fb36a803ff1e94a6 100644
--- a/directional.cpp
+++ b/directional.cpp
@@ -58,7 +58,7 @@ int directional(const GmshFem& gmshFem)
ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
- ConfigurationInterface* configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
unsigned int n_freq_tot = 0;
if(!gmshFem.userDefinedParameter(n_freq_tot, "n_freq"))
@@ -85,7 +85,7 @@ int directional(const GmshFem& gmshFem)
msg::warning << "filename #"+suffix+" = "<< filename[f] << msg::endl;
}
}
- const Data<T_Physic> d0(filename,frequency,configuration);
+ const Data<T_Physic> d0(filename,frequency,configuration.get());
unsigned int n_group = 0;
if(!gmshFem.userDefinedParameter(n_group, "n_group"))
@@ -139,7 +139,7 @@ int directional(const GmshFem& gmshFem)
}
const Data<T_Physic> d0g(freq_idx,d0);
- ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(d0g,configuration,gmshFem);
+ auto objective = newObjective<T_Physic>(d0g,configuration.get(),gmshFem);
std::vector<EquationInterface<T_Physic>*> pequation(n_freq);
std::vector<double> subfrequency(n_freq);
@@ -149,10 +149,10 @@ int directional(const GmshFem& gmshFem)
subfrequency[f] = frequency[freq_idx[f]];
wave::Discretization<T_Physic> w_discret(gmshFem,suffix_f);
- pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*subfrequency[f],configuration,w_discret,gmshFem,suffix_f);
+ pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*subfrequency[f],configuration.get(),w_discret,gmshFem,suffix_f);
}
- StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,nullptr,nullptr,subfrequency,pequation,objective);
+ StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration.get(),nullptr,nullptr,subfrequency,pequation,objective.get());
FunctionalInterface* const functional = statefunctional;
double eps;
@@ -193,7 +193,7 @@ int directional(const GmshFem& gmshFem)
}
ScalarPiecewiseFunction< std::complex< double > > spf_dm;
- spf_dm.addFunction(Redm+im*Imdm,configuration->model_unknown(Support::BLK) | configuration->model_unknown(Support::BND));
+ spf_dm.addFunction(Redm+im*Imdm,configuration.get()->model_unknown(Support::BLK) | configuration->model_unknown(Support::BND));
dm[c] = spf_dm;
@@ -344,11 +344,9 @@ int directional(const GmshFem& gmshFem)
{
delete pequation[f];
}
- delete objective;
msg::unindent();
}
- delete configuration;
delete parametrization;
return 0;
diff --git a/directional_field.cpp b/directional_field.cpp
index 20d08f5dc78d2b58a61a000cfe34102550ae5ea4..b52403046f978fa7277e84796ade8947e4c7f622 100644
--- a/directional_field.cpp
+++ b/directional_field.cpp
@@ -32,7 +32,7 @@ int directional(const GmshFem& gmshFem)
ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
- ConfigurationInterface* configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
unsigned int n_freq_tot = 0;
if(!gmshFem.userDefinedParameter(n_freq_tot, "n_freq"))
@@ -59,7 +59,7 @@ int directional(const GmshFem& gmshFem)
msg::warning << "filename #"+suffix+" = "<< filename[f] << msg::endl;
}
}
- const Data<T_Physic> d0(filename,frequency,configuration);
+ const Data<T_Physic> d0(filename,frequency,configuration.get());
unsigned int n_group = 0;
if(!gmshFem.userDefinedParameter(n_group, "n_group"))
@@ -113,7 +113,7 @@ int directional(const GmshFem& gmshFem)
}
const Data<T_Physic> d0g(freq_idx,d0);
- ObjectiveInterface<T_Physic>* const objective = (d0g.nf()!=0 ? newObjective<T_Physic>(d0g,configuration,gmshFem) : nullptr);
+ std::unique_ptr<ObjectiveInterface<T_Physic>> objective = (d0g.nf()!=0 ? newObjective<T_Physic>(d0g,configuration.get(),gmshFem) : nullptr);
std::vector<EquationInterface<T_Physic>*> pequation(n_freq);
std::vector<double> subfrequency(n_freq);
@@ -123,16 +123,16 @@ int directional(const GmshFem& gmshFem)
subfrequency[f] = frequency[freq_idx[f]];
wave::Discretization<T_Physic> w_discret(gmshFem,suffix_f);
- pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*subfrequency[f],configuration,w_discret,gmshFem,suffix_f);
+ pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*subfrequency[f],configuration.get(),w_discret,gmshFem,suffix_f);
}
model::Discretization m_discret(gmshFem,suffix_g);
- ModelFieldSweeperInterface* sweeper_m = newModelFieldSweeper("m",configuration,m_discret,gmshFem);
- ModelFieldSweeperInterface* sweeper_dm = newModelFieldSweeper("dm",configuration,m_discret,gmshFem);
+ ModelFieldSweeperInterface* sweeper_m = newModelFieldSweeper("m",configuration.get(),m_discret,gmshFem);
+ ModelFieldSweeperInterface* sweeper_dm = newModelFieldSweeper("dm",configuration.get(),m_discret,gmshFem);
- RegularizationInterface* const regularization = newRegularization(configuration,(*sweeper_m)(0,0), gmshFem);
+ RegularizationInterface* const regularization = newRegularization(configuration.get(),(*sweeper_m)(0,0), gmshFem);
- StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,nullptr,regularization,subfrequency,pequation,objective);
+ StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration.get(),nullptr,regularization,subfrequency,pequation,objective.get());
FunctionalInterface* const functional = statefunctional;
@@ -214,12 +214,10 @@ int directional(const GmshFem& gmshFem)
{
delete pequation[f];
}
- delete objective;
delete regularization;
msg::unindent();
}
- delete configuration;
delete parametrization;
return 0;
diff --git a/error.cpp b/error.cpp
index c4896285ffbb29f60eebcf8478413957d63c1534..d37cfb736143656580fa2405eb79dabbe43e43de 100644
--- a/error.cpp
+++ b/error.cpp
@@ -42,7 +42,7 @@ int error(const GmshFem& gmshFem)
const ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
- ConfigurationInterface* configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
unsigned int write_error_field = 0;
gmshFem.userDefinedParameter(write_error_field, "write_error_field");
diff --git a/gradient.cpp b/gradient.cpp
index bb6bb2b33dcc895070360d3e4cb0fdbdf4792bc2..b71b1819ed4165b8e1ec57b6a8b4ecb20b91ebb3 100644
--- a/gradient.cpp
+++ b/gradient.cpp
@@ -45,7 +45,7 @@ int gradient(const GmshFem& gmshFem)
const ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
- const ConfigurationInterface* const configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
unsigned int n_freq_tot = 0;
@@ -73,7 +73,7 @@ int gradient(const GmshFem& gmshFem)
msg::warning << "filename #"+suffix+" = "<< filename[f] << msg::endl;
}
}
- const Data<T_Physic> d0(filename,frequency,configuration);
+ const Data<T_Physic> d0(filename,frequency,configuration.get());
for (unsigned int c = 0; c < parametrization->size(); c++)
{
@@ -81,7 +81,7 @@ int gradient(const GmshFem& gmshFem)
gmshfem::post::save(configuration->m0()[c], configuration->model_unknown(Support::BND), name+"_m_bndc"+std::to_string(c)+"f0", "pos", "");
}
model::Discretization m_discret(gmshFem);
- sobolev_l2::InnerProduct innerproduct(configuration,m_discret,gmshFem);
+ sobolev_l2::InnerProduct innerproduct(configuration.get(),m_discret,gmshFem);
unsigned int n_group = 0;
if(!gmshFem.userDefinedParameter(n_group, "n_group"))
@@ -134,7 +134,7 @@ int gradient(const GmshFem& gmshFem)
}
const Data<T_Physic> d0g(freq_idx,d0);
- ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(d0g,configuration,gmshFem);
+ auto objective = newObjective<T_Physic>(d0g,configuration.get(),gmshFem);
std::vector<EquationInterface<T_Physic>*> pequation(n_freq);
std::vector<double> subfrequency(n_freq);
@@ -144,10 +144,10 @@ int gradient(const GmshFem& gmshFem)
subfrequency[f] = frequency[freq_idx[f]];
wave::Discretization<T_Physic> w_discret(gmshFem,suffix_f);
- pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*subfrequency[f],configuration,w_discret,gmshFem,suffix_f);
+ pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*subfrequency[f],configuration.get(),w_discret,gmshFem,suffix_f);
}
- StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,&innerproduct,nullptr,subfrequency,pequation,objective);
+ StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration.get(),&innerproduct,nullptr,subfrequency,pequation,objective.get());
FunctionalInterface* const functional = statefunctional;
functional->setModel(configuration->m0());
@@ -203,14 +203,12 @@ int gradient(const GmshFem& gmshFem)
msg::unindent();
delete statefunctional;
- delete objective;
for (unsigned int f = 0; f < n_freq; f++)
{
delete pequation[f];
}
msg::unindent();
}
- delete configuration;
delete parametrization;
return 0;
diff --git a/inversion.cpp b/inversion.cpp
index 761006612c8a4ad56520d6755db82d3b730a09f5..af1cfab517475e16c3a2945bcf918bb1ba5ad8fd 100644
--- a/inversion.cpp
+++ b/inversion.cpp
@@ -28,7 +28,7 @@ int inversion(const GmshFem& gmshFem)
const ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
- const ConfigurationInterface* const configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
for (unsigned int c = 0; c < parametrization->size(); c++)
{
@@ -61,7 +61,7 @@ int inversion(const GmshFem& gmshFem)
msg::warning << "filename #"+suffix+" = "<< filename[f] << msg::endl;
}
}
- const Data<T_Physic> d0(filename,frequency,configuration);
+ const Data<T_Physic> d0(filename,frequency,configuration.get());
unsigned int n_group = 0;
if(!gmshFem.userDefinedParameter(n_group, "n_group"))
@@ -132,17 +132,17 @@ int inversion(const GmshFem& gmshFem)
double m0_scale=0.;
if(gmshFem.userDefinedParameter(m0_scale, "m0_scale"))
{
- m = laplace_filter(configuration->m0(),m0_scale,configuration, m_discret, integrationType);
+ m = laplace_filter(configuration->m0(),m0_scale,configuration.get(), m_discret, integrationType);
}
else
{
- m = rediscretize(configuration->m0(), configuration, m_discret, integrationType);
+ m = rediscretize(configuration->m0(), configuration.get(), m_discret, integrationType);
}
m.write(name+"_m_initial", "pos", "");
}
else
{
- m = rediscretize(m,configuration, m_discret, integrationType);
+ m = rediscretize(m,configuration.get(), m_discret, integrationType);
}
const Data<T_Physic> d0g(freq_idx,d0);
@@ -152,10 +152,10 @@ int inversion(const GmshFem& gmshFem)
std::string suffix_f = std::to_string(freq_idx[f]);
wave::Discretization<T_Physic> w_discret(gmshFem,suffix_f);
- pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*frequency[freq_idx[f]],configuration,w_discret,gmshFem,suffix_f);
+ pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*frequency[freq_idx[f]],configuration.get(),w_discret,gmshFem,suffix_f);
}
- GlobalMinimumSearchInterface<T_Physic>* const globalminimumsearch = newGlobalMinimumSearch<T_Physic>(name,configuration, gmshFem, suffix_g);
+ GlobalMinimumSearchInterface<T_Physic>* const globalminimumsearch = newGlobalMinimumSearch<T_Physic>(name,configuration.get(), gmshFem, suffix_g);
(*globalminimumsearch)(&m,d0g,pequation);
m.write(name+"_mg"+suffix_g, "pos", "");
@@ -170,7 +170,6 @@ int inversion(const GmshFem& gmshFem)
msg::unindent();
}
- delete configuration;
delete parametrization;
return 0;
}
diff --git a/ip_comparison.cpp b/ip_comparison.cpp
index 969785a6b7cf649b1ec37b5ace3f3e880d673c24..b717c27b44c9a9cf76356a1b7e03408be2aec981 100644
--- a/ip_comparison.cpp
+++ b/ip_comparison.cpp
@@ -32,7 +32,7 @@ int ip_comparison(const GmshFem& gmshFem)
const ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
- const ConfigurationInterface* const configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
std::string integrationType;
gmshFem.userDefinedParameter(integrationType, "integration_type");
@@ -67,7 +67,7 @@ int ip_comparison(const GmshFem& gmshFem)
msg::warning << "filename #"+suffix+" = "<< filename[f] << msg::endl;
}
}
- const Data<T_Physic> d0(filename,frequency,configuration);
+ const Data<T_Physic> d0(filename,frequency,configuration.get());
unsigned int n_group = 0;
if(!gmshFem.userDefinedParameter(n_group, "n_group"))
@@ -123,10 +123,10 @@ int ip_comparison(const GmshFem& gmshFem)
msg::print << "Perform comparison" << msg::endl;
msg::indent();
- ModelField m = laplace_filter(configuration->m0(),filter_scale,configuration,m_discret,integrationType);
+ ModelField m = laplace_filter(configuration.get()->m0(),filter_scale,configuration.get(),m_discret,integrationType);
m.write(name+"_m_smooth"+suffix_g);
- ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(d0g,configuration,gmshFem);
+ auto objective = newObjective<T_Physic>(d0g,configuration.get(),gmshFem);
std::vector<EquationInterface<T_Physic>*> pequation(n_freq);
std::vector<double> subfrequency(n_freq);
@@ -136,12 +136,12 @@ int ip_comparison(const GmshFem& gmshFem)
subfrequency[f] = frequency[freq_idx[f]];
wave::Discretization<T_Physic> w_discret(gmshFem,suffix_f);
- pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*subfrequency[f],configuration,w_discret,gmshFem,suffix_f);
+ pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*subfrequency[f],configuration.get(),w_discret,gmshFem,suffix_f);
}
- InnerProductInterface* const innerproduct = newInnerProduct(configuration,m_discret,gmshFem,suffix_g);
- RegularizationInterface* const regularization = newRegularization(configuration,m,gmshFem,suffix_g);
- StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,innerproduct,regularization,subfrequency,pequation,objective);
+ auto innerproduct = newInnerProduct(configuration.get(),m_discret,gmshFem,suffix_g);
+ RegularizationInterface* const regularization = newRegularization(configuration.get(),m,gmshFem,suffix_g);
+ StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration.get(),innerproduct.get(),regularization,subfrequency,pequation,objective.get());
statefunctional->setModel(m);
double j = statefunctional->performance();
msg::print << "performance = " << j << msg::endl;
@@ -206,17 +206,14 @@ int ip_comparison(const GmshFem& gmshFem)
msg::unindent();
delete statefunctional;
- delete objective;
for (unsigned int f = 0; f < n_freq; f++)
{
delete pequation[f];
}
- delete innerproduct;
delete regularization;
msg::unindent();
}
msg::unindent();
- delete configuration;
delete parametrization;
return 0;
diff --git a/multiscale.cpp b/multiscale.cpp
index 8982bac61a5c3acc734c05ca92067514e1d27a37..a8b335e0ff048a94ab1feb5f506bb18ab77f795d 100644
--- a/multiscale.cpp
+++ b/multiscale.cpp
@@ -33,7 +33,7 @@ int multiscale(const GmshFem& gmshFem)
name += suffix;
ParametrizationInterface* const parametrization = newParametrization<Physic::acoustic>(gmshFem);
- const ConfigurationInterface* const configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
std::string integrationType;
gmshFem.userDefinedParameter(integrationType, "model_IntegrationType");
@@ -68,9 +68,9 @@ int multiscale(const GmshFem& gmshFem)
msg::warning << "filename #"+suffix+" = "<< filename[f] << msg::endl;
}
}
- const Data<T_Physic> d0(filename,frequency,configuration);
+ const Data<T_Physic> d0(filename,frequency,configuration.get());
- ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(d0,configuration,gmshFem);
+ auto objective = newObjective<T_Physic>(d0,configuration.get(),gmshFem);
std::vector<EquationInterface<T_Physic>*> pequation(n_freq);
for (unsigned int f = 0; f < n_freq; f++)
@@ -78,9 +78,9 @@ int multiscale(const GmshFem& gmshFem)
std::string suffix_f = std::to_string(f);
wave::Discretization<T_Physic> w_discret(gmshFem,suffix_f);
- pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*frequency[f],configuration,w_discret,gmshFem,suffix_f);
+ pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,f,2.*M_PI*frequency[f],configuration.get(),w_discret,gmshFem,suffix_f);
}
- StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,nullptr,nullptr,frequency,pequation,objective);
+ StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration.get(),nullptr,nullptr,frequency,pequation,objective.get());
FunctionalInterface* const functional = statefunctional;
msg::print << "Perform multiscale analysis" << msg::endl;
@@ -124,13 +124,11 @@ int multiscale(const GmshFem& gmshFem)
delete statefunctional;
- delete objective;
for (unsigned int f = 0; f < n_freq; f++)
{
delete pequation[f];
}
delete sweeper;
- delete configuration;
delete parametrization;
return 0;
diff --git a/ob_comparison.cpp b/ob_comparison.cpp
index daa8439cdd60d415702e64eedc0f5e3318feb1a2..841fa2b401716071e179674e1cb9efe676a0de7d 100644
--- a/ob_comparison.cpp
+++ b/ob_comparison.cpp
@@ -41,7 +41,7 @@ int ob_comparison(const GmshFem& gmshFem)
ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
- ConfigurationInterface* configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
std::string integrationType;
gmshFem.userDefinedParameter(integrationType, "model_IntegrationType");
@@ -51,7 +51,7 @@ int ob_comparison(const GmshFem& gmshFem)
double filter_scale = 0.;
gmshFem.userDefinedParameter(filter_scale, "filter_scale");
- ModelUpdater mu(configuration,nullptr,nullptr);
+ ModelUpdater mu(configuration.get(),nullptr,nullptr);
unsigned int n_freq = 1;
gmshFem.userDefinedParameter(n_freq, "n_freq");
@@ -69,22 +69,22 @@ int ob_comparison(const GmshFem& gmshFem)
msg::print << "--- Frequency #"+suffix+" = " << frequency << " --- " << msg::endl;
msg::indent();
- ModelField m = laplace_filter(configuration->m0(),filter_scale,configuration,m_discret,integrationType);
+ ModelField m = laplace_filter(configuration->m0(),filter_scale,configuration.get(),m_discret,integrationType);
m.write(name+"_m_smooth"+suffix);
wave::Discretization<T_Physic> w_discret(gmshFem,suffix);
parametrization->pulsation(2.*M_PI*frequency);
- EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem);
- EquationInterface<T_Physic>* const pequation = new ParametrizedEquation<T_Physic>(configuration,equation,parametrization,gmshFem);
+ EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration.get(),w_discret,gmshFem);
+ EquationInterface<T_Physic>* const pequation = new ParametrizedEquation<T_Physic>(configuration.get(),equation,parametrization,gmshFem);
std::string filename = "noname";
gmshFem.userDefinedParameter(filename, "data"+suffix);
- Data<T_Physic> d0(filename,configuration);
+ Data<T_Physic> d0(filename,configuration.get());
- continuousl2distance::Objective<T_Physic>* const objective = new continuousl2distance::Objective<T_Physic>(&d0,configuration,gmshFem,suffix);
+ continuousl2distance::Objective<T_Physic>* const objective = new continuousl2distance::Objective<T_Physic>(&d0,configuration.get(),gmshFem,suffix);
- WaveUpdater<T_Physic> wu(configuration,nullptr,pequation);
- DataUpdater<T_Physic> du(configuration,&wu,objective);
+ WaveUpdater<T_Physic> wu(configuration.get(),nullptr,pequation);
+ DataUpdater<T_Physic> du(configuration.get(),&wu,objective);
std::array<bool,5> dataNeedUpToDate = {true,false,false,false,false};
std::array<bool,4> modelNeedUpToDate = {true,false,false,false};
Data<T_Physic> d = (du.get(dataNeedUpToDate,mu.get(modelNeedUpToDate,&m))).state(Type::FS);
@@ -150,7 +150,6 @@ int ob_comparison(const GmshFem& gmshFem)
msg::unindent();
}
msg::unindent();
- delete configuration;
delete parametrization;
return 0;
diff --git a/preconditioner.cpp b/preconditioner.cpp
index 04dd0623cc14ff935d64e921290a6f4c41ff6997..83e88cc0189c7a5e574e15c23f29889cc0ade888 100644
--- a/preconditioner.cpp
+++ b/preconditioner.cpp
@@ -27,7 +27,7 @@ int preconditioner(const GmshFem& gmshFem)
const ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
- const ConfigurationInterface* const configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
unsigned int n_freq_tot = 0;
if(!gmshFem.userDefinedParameter(n_freq_tot, "n_freq"))
@@ -35,7 +35,7 @@ int preconditioner(const GmshFem& gmshFem)
throw Exception("Total frequency number could not be found.");
}
std::vector<double> frequency(n_freq_tot);
- const Data<T_Physic>* const d0 = new Data<T_Physic>(configuration);
+ const Data<T_Physic>* const d0 = new Data<T_Physic>(configuration.get());
for (unsigned int f = 0; f < n_freq_tot; f++)
{
std::string suffix = std::to_string(f);
@@ -83,7 +83,7 @@ int preconditioner(const GmshFem& gmshFem)
msg::print << "Compute preconditioner" << msg::endl;
- ModelUpdater mu(configuration,nullptr,nullptr);
+ ModelUpdater mu(configuration.get(),nullptr,nullptr);
ModelFunction m0 = configuration->m0();
std::array<bool,4> needUpToDate = {true,false,false,false};
mu.get(needUpToDate,&m0).write(Type::FS,name+"_m"+suffix_g);
@@ -100,14 +100,14 @@ int preconditioner(const GmshFem& gmshFem)
std::string suffix_f = std::to_string(freq_idx[f]);
wave::Discretization<T_Physic> w_discret(gmshFem,suffix_f);
- equation[f] = newEquation(2.*M_PI*frequency[freq_idx[f]],configuration,w_discret,gmshFem,suffix_f);
- pequation[f] = new ParametrizedEquation<T_Physic>(configuration,equation[f],parametrization,gmshFem,suffix_f);
+ equation[f] = newEquation(2.*M_PI*frequency[freq_idx[f]],configuration.get(),w_discret,gmshFem,suffix_f);
+ pequation[f] = new ParametrizedEquation<T_Physic>(configuration.get(),equation[f],parametrization,gmshFem,suffix_f);
- objective[f] = newObjective<T_Physic>(d0,configuration,gmshFem,suffix_f);
+ objective[f] = newObjective<T_Physic>(d0,configuration.get(),gmshFem,suffix_f);
objective[f]->link(&mu,&du[f]);
- new(&wu[f]) WaveUpdater<T_Physic>(configuration,nullptr,pequation[f]);
- new(&du[f]) DataUpdater<T_Physic>(configuration,&wu[f],objective[f]);
+ new(&wu[f]) WaveUpdater<T_Physic>(configuration.get(),nullptr,pequation[f]);
+ new(&du[f]) DataUpdater<T_Physic>(configuration.get(),&wu[f],objective[f]);
sue[f] = new SensitivityUpdaterFromEquation<T_Physic>(&wu[f],&du[f],pequation[f]);
}
SensitivityUpdaterFromVector su(n_freq,sue);
@@ -130,7 +130,6 @@ int preconditioner(const GmshFem& gmshFem)
msg::unindent();
}
msg::unindent();
- delete configuration;
delete parametrization;
return 0;
diff --git a/specific/configuration/circular_acquisition.cpp b/specific/configuration/circular_acquisition.cpp
index 756729138122bc85c489d051310c7d8bf09ab061..59002d6e213c44bef28166f4a77f4197bf99881b 100644
--- a/specific/configuration/circular_acquisition.cpp
+++ b/specific/configuration/circular_acquisition.cpp
@@ -267,26 +267,7 @@ namespace circular_acquisition
}
}
- void Configuration::mesh() const
- {
- msg::print << "Generate meshes" << msg::endl;
- gmsh::option::setNumber("General.Terminal", 1);
-
- gmodel::add(_name);
- if(!_remesh)
- {
- gmsh::open(_name+".msh");
- return;
- }
- wave_mesh();
- data_mesh();
-
- factory::synchronize();
- gmodel::mesh::generate();
-
- gmsh::write(_name+".msh");
- }
void Configuration::wave_mesh() const
{
int p0 = factory::addPoint(0., 0., 0., _h);
diff --git a/specific/configuration/circular_acquisition.h b/specific/configuration/circular_acquisition.h
index 40842cf1cd5da07e07c49c5feffc0f914b74bbfb..75e4f3de4eead27fbaa761a62a12fa4c5616c495 100644
--- a/specific/configuration/circular_acquisition.h
+++ b/specific/configuration/circular_acquisition.h
@@ -70,7 +70,6 @@ namespace circular_acquisition
bool areFilled() const {return _areFilled;};
unsigned int ner() const {return _np;};
- virtual void mesh() const;
virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const;
virtual std::array<double,2> index_to_data_coordinate(unsigned int s, unsigned int r) const;
diff --git a/specific/configuration/flexible_acquisition.cpp b/specific/configuration/flexible_acquisition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3850bb1d1eb83971fbfcc1febee4c7c5faeb2722
--- /dev/null
+++ b/specific/configuration/flexible_acquisition.cpp
@@ -0,0 +1,457 @@
+//Standard Library
+#include <numeric>
+#include <array>
+
+//GmshFEM Library
+#include "GmshFem.h"
+#include "Exception.h"
+#include "Message.h"
+#include "Domain.h"
+#include "Function.h"
+#include "CSVio.h"
+
+//Gmsh Library
+#include "gmsh.h"
+
+//GmshFWI Library
+#include "../../common/model/macro.h"
+#include "../../common/data/element.h"
+#include "../wave/correlation.h"
+#include "flexible_acquisition.h"
+#include "yaml_interface.h"
+
+namespace gmodel = gmsh::model;
+namespace factory = gmsh::model::geo;
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::domain;
+using namespace gmshfem::function;
+
+static const std::complex< double > im = std::complex< double >(0., 1.);
+
+namespace flexible_acquisition
+{
+
+ /*
+ * class Configuration
+ */
+ Configuration::Configuration(std::string name, const ParametrizationInterface* const parametrization, const GmshFem& gmshFem) : ConfigurationInterface(name, parametrization, gmshFem)
+ {
+ msg::print << "Initialize flexible configuration" << msg::endl;
+
+ /*
+ * MESH
+ */
+ if(
+ !(
+ gmshFem.userDefinedParameter(_H, "H") &&
+ gmshFem.userDefinedParameter(_L, "L") &&
+ gmshFem.userDefinedParameter(_Hsub, "Hsub") &&
+ gmshFem.userDefinedParameter(_h, "h")
+ )
+ )
+ {
+ throw common::Exception("A geometric parameter could not be found.");
+ }
+
+ // SETUP FROM YAML
+ std::string shotsPath;
+ if (!gmshFem.userDefinedParameter(shotsPath, "shots_config"))
+ {
+ throw common::Exception("Path to a shots_config is needed in flexible_acquisition.");
+ }
+ ShotsConfigurationYAML shotsConfig(shotsPath);
+ _ns = shotsConfig.numShots();
+
+ // Configure points
+ for (const auto &coords : shotsConfig.points())
+ {
+ _er_positions.push_back(coords);
+ }
+
+ _np = _er_positions.size();
+
+
+ mesh();
+
+
+ // Configure shots
+ for (unsigned s = 0; s < shotsConfig.numShots(); ++s)
+ {
+ _emitter.push_back(shotsConfig.emitters(s));
+ _receiver.push_back(shotsConfig.receivers(s));
+ }
+
+ for (unsigned p = 0; p < _er_positions.size(); ++p) {
+ _point.push_back("emitter_receiver_"+std::to_string(p));
+ _points |= _point[p];
+ }
+
+
+
+ /*
+ * DOMAIN
+ */
+ // TODO: actual data mesh!
+ _data_omega = Domain("supersurface") | Domain("subsurface");
+
+ _subsurface[Support::BLK] = Domain("subsurface");
+ _subsurface[Support::BND] = Domain("subsurface_bnd");
+ _supersurface[Support::BLK] = Domain("supersurface");
+ _supersurface[Support::BND] = Domain("supersurface_bnd");
+
+ std::string unknownRegion = "subsurface";
+ gmshFem.userDefinedParameter(unknownRegion, "unknown");
+ if (unknownRegion == "subsurface")
+ {
+ _model_known[Support::BLK] = _supersurface[Support::BLK];
+ _model_known[Support::BND] = _supersurface[Support::BND];
+
+ _model_unknown[Support::BLK] = _subsurface[Support::BLK];
+ _model_unknown[Support::BND] = _subsurface[Support::BND];
+ }
+ else if (unknownRegion == "none")
+ {
+ _model_known[Support::BLK] = _supersurface[Support::BLK] | _subsurface[Support::BLK];
+ _model_known[Support::BND] = _supersurface[Support::BND] | _subsurface[Support::BND];
+ }
+ else if (unknownRegion == "all")
+ {
+ _model_unknown[Support::BLK] = _supersurface[Support::BLK] | _subsurface[Support::BLK];
+ _model_unknown[Support::BND] = _supersurface[Support::BND] | _subsurface[Support::BND];
+ }
+ else {
+ throw Exception("Invalid unknown region type: " + unknownRegion);
+ }
+
+ _wave_omega[Support::BLK] = _model_known[Support::BLK] | _model_unknown[Support::BLK];
+ _wave_omega[Support::BND] = _model_known[Support::BND] | _model_unknown[Support::BND];
+
+ // Once Dirichlet BCs are added
+ //_named_domains() ...
+ _named_domains["top_bnd"] = Domain("top_bnd");
+
+ /*
+ * Reference model function
+ */
+ // TODO: remove debugging stuff
+ setupInitialModel(parametrization, gmshFem);
+ gmshfem::post::save(_m0[0], _model_known[Support::BLK], "model_known");
+ gmshfem::post::save(_m0[0], _model_unknown[Support::BLK], "model_unknown");
+
+
+ }
+
+
+ void Configuration::wave_mesh() const
+ {
+ int p0 = factory::addPoint(0., -_H, 0., _h);
+ int p1 = factory::addPoint(_L, -_H, 0., _h);
+ int p2 = factory::addPoint(_L, -_Hsub, 0., _h);
+ int p3 = factory::addPoint(0., -_Hsub, 0., _h);
+ int p4 = factory::addPoint(_L, 0., 0., _h);
+ int p5 = factory::addPoint(0., 0., 0., _h);
+
+ // All horizontal lines defined left to right, vertical bottom to top
+ int lineBottom = factory::addLine(p0, p1);
+ int lineSub = factory::addLine(p3, p2);
+ int lineTop = factory::addLine(p5, p4);
+ int lineLeftSub = factory::addLine(p0, p3);
+ int lineLeftSuper = factory::addLine(p3, p5);
+ int lineRightSub = factory::addLine(p1, p2);
+ int lineRightSuper = factory::addLine(p2, p4);
+
+ int clSub = factory::addCurveLoop({lineBottom, lineRightSub, -lineSub, -lineLeftSub});
+ int clSuper = factory::addCurveLoop({lineSub, lineRightSuper, -lineTop, -lineLeftSuper});
+
+ int sSub = factory::addPlaneSurface({clSub});
+ int sSuper = factory::addPlaneSurface({clSuper});
+
+ std::vector<int> p_er;
+ // Emitters and receivers
+ for (auto pos: _er_positions) {
+ p_er.push_back(factory::addPoint(pos.first, pos.second, 0.0, _h));
+ }
+
+ factory::synchronize();
+
+ // Embed the e/r. Check if they're in the subsurface or supersurface
+ for (unsigned i = 0; i < p_er.size(); ++i) {
+ bool isSub = _er_positions[i].second < (-_Hsub);
+ gmodel::mesh::embed(0, {p_er[i]}, 2, isSub ? sSub : sSuper);
+
+ }
+
+ gmodel::addPhysicalGroup(2, {sSub}, 21);
+ gmodel::setPhysicalName(2, 21, "subsurface");
+ gmodel::addPhysicalGroup(2, {sSuper}, 22);
+ gmodel::setPhysicalName(2, 22, "supersurface");
+
+ gmodel::addPhysicalGroup(1, {lineLeftSub, lineBottom, lineRightSub}, 11);
+ gmodel::setPhysicalName(1, 11, "subsurface_bnd");
+ gmodel::addPhysicalGroup(1, {lineLeftSuper, lineTop, lineRightSuper}, 12);
+ gmodel::setPhysicalName(1, 12, "supersurface_bnd");
+
+ gmodel::addPhysicalGroup(1, {lineTop}, 13);
+ gmodel::setPhysicalName(1, 13, "top_bnd");
+
+ for (unsigned p = 0; p < p_er.size(); ++p) {
+ gmodel::addPhysicalGroup(0, {p_er[p]}, 1000+p);
+ gmodel::setPhysicalName(0, 1000+p, "emitter_receiver_"+std::to_string(p));
+ }
+
+ }
+ void Configuration::data_mesh() const
+ {
+ // TODO
+ /*
+ double DAngle = 2. * M_PI /((double) _np);
+ double Ds_e = DAngle * _emitter_skip * _rer;
+ double Ds_r = DAngle * _receiver_skip * _rer;
+ double Le0 = DAngle * _emitter_offset * _rer;
+ double Lr0 = DAngle * _receiver_offset * _rer;
+ double Le = (_ne-1) * Ds_e;
+ double Lr = (_nr-1) * Ds_r;
+
+ int p0 = factory::addPoint(Le0, Lr0, 1.);
+ int p1 = factory::addPoint(Le0+Le, Lr0, 1.);
+ int p2 = factory::addPoint(Le0+Le, Lr0+Lr, 1.);
+ int p3 = factory::addPoint(Le0, Lr0+Lr, 1.);
+
+ int l0 = factory::addLine(p0, p1);
+ int l1 = factory::addLine(p1, p2);
+ int l2 = factory::addLine(p2, p3);
+ int l3 = factory::addLine(p3, p0);
+
+ int cl = factory::addCurveLoop({l0,l1,l2,l3});
+ int s1 = factory::addPlaneSurface({cl});
+
+ factory::mesh::setTransfiniteCurve(l0,_ne);
+ factory::mesh::setTransfiniteCurve(l1,_nr);
+ factory::mesh::setTransfiniteCurve(l2,_ne);
+ factory::mesh::setTransfiniteCurve(l3,_nr);
+ factory::mesh::setTransfiniteSurface(s1,"Left",{p0,p1,p2,p3});
+
+ factory::mesh::setRecombine(2,s1);
+
+ gmodel::addPhysicalGroup(2, {s1}, 1);
+ gmodel::setPhysicalName(2, 1, "data_omega");
+
+ */
+ }
+
+ void Configuration::setupInitialModel(const ParametrizationInterface* const parametrization, const GmshFem& gmshFem) {
+ /*
+ * Reference model function
+ */
+ _m_super.resize(model_size());
+ _m_sub.resize(model_size());
+ for (unsigned int c = 0; c < model_size(); c++)
+ {
+ ScalarPiecewiseFunction<std::complex<double>> m0;
+ std::string suffix = "c" + std::to_string(c);
+ double Rem_super, Imm_super, Rem_sub, Imm_sub;
+ if (!(
+ gmshFem.userDefinedParameter(Rem_super, "Re(m_super" + suffix + ")") && gmshFem.userDefinedParameter(Imm_super, "Im(m_super" + suffix + ")") &&
+ gmshFem.userDefinedParameter(Rem_sub, "Re(m_sub" + suffix + ")") && gmshFem.userDefinedParameter(Imm_sub, "Im(m_sub" + suffix + ")")))
+ {
+ throw common::Exception("A model component " + suffix + " could not be found.");
+ }
+ _m_super[c] = Rem_super + im * Imm_super;
+ _m_sub[c] = Rem_sub + im * Imm_sub;
+ _mc[c] = _m_sub[c];
+
+ m0.addFunction(_m_super[c], _supersurface[Support::BLK] | _supersurface[Support::BND] | _points);
+
+ std::string m0_type;
+ if (!gmshFem.userDefinedParameter(m0_type, "m0_type" + suffix))
+ {
+ throw Exception("Reference model type could not be found.");
+ }
+
+ if (m0_type == "file")
+ {
+ std::string path = "";
+ if (!gmshFem.userDefinedParameter(path, "m0_path" + suffix))
+ {
+ throw common::Exception("Path to subsurface data could not be found.");
+ }
+ m0.addFunction(bilinearInterpolation(path), _subsurface[Support::BLK] | _subsurface[Support::BND]);
+ }
+ else if (m0_type == "file.pos")
+ {
+ std::string path = "";
+ if (!gmshFem.userDefinedParameter(path, "m0_path" + suffix))
+ {
+ throw common::Exception("Path to subsurface data could not be found.");
+ }
+
+ gmsh::merge(path + suffix + ".pos");
+ // Fix for weird tags when loading .pos files.
+ std::vector<int> tags;
+ gmsh::view::getTags(tags);
+ auto view = tags.back();
+ ScalarFunction<std::complex<double>> mpos = probeScalarView<std::complex<double>>(view);
+
+ m0.addFunction(mpos, _subsurface[Support::BLK] | _subsurface[Support::BND]);
+ }
+ /*
+ else if (m0_type == "inverse_linear_squared")
+ {
+ double Rea_0, Ima_0;
+ if (!(
+ gmshFem.userDefinedParameter(Rea_0, "Re(a0" + suffix + ")") && gmshFem.userDefinedParameter(Ima_0, "Im(a0" + suffix + ")")))
+ {
+ throw common::Exception("Initial model parameter (a0) could not be found.");
+ }
+
+ double Rea_H, Ima_H;
+ if (!(
+ gmshFem.userDefinedParameter(Rea_H, "Re(aH" + suffix + ")") && gmshFem.userDefinedParameter(Ima_H, "Im(aH" + suffix + ")")))
+ {
+ throw common::Exception("Initial model parameter (aH) could not be found.");
+ }
+
+ ScalarFunction<std::complex<double>> num = (Rea_0 + im * Ima_0) - ((Rea_H + im * Ima_H) - (Rea_0 + im * Ima_0)) / H() * y<std::complex<double>>();
+ m0.addFunction(1. / pow(num, 2), _subsurface[Support::BLK] | _subsurface[Support::BND]);
+ }
+ else if (m0_type == "linear")
+ {
+ double Rea_0, Ima_0;
+ if (!(
+ gmshFem.userDefinedParameter(Rea_0, "Re(a0" + suffix + ")") && gmshFem.userDefinedParameter(Ima_0, "Im(a0" + suffix + ")")))
+ {
+ throw common::Exception("Initial model parameter (a0) could not be found.");
+ }
+
+ double Rea_H, Ima_H;
+ if (!(
+ gmshFem.userDefinedParameter(Rea_H, "Re(aH" + suffix + ")") && gmshFem.userDefinedParameter(Ima_H, "Im(aH" + suffix + ")")))
+ {
+ throw common::Exception("Initial model parameter (aH) could not be found.");
+ }
+
+ ScalarFunction<std::complex<double>> lin = (Rea_0 + im * Ima_0) - ((Rea_H + im * Ima_H) - (Rea_0 + im * Ima_0)) / H() * (y<std::complex<double>>() + _ymin);
+ m0.addFunction(lin, _subsurface[Support::BLK] | _subsurface[Support::BND]);
+ }*/
+ else if (m0_type == "constant")
+ {
+ m0.addFunction(_m_sub[c], _subsurface[Support::BLK] | _subsurface[Support::BND]);
+ }
+ else
+ {
+ throw common::Exception("Initial model type (" + suffix + ") " + m0_type + " is unknown.");
+ }
+ _m0.push_back(m0);
+ }
+ }
+
+ std::array<unsigned int,2> Configuration::data_coordinate_to_index(double xs, double xr) const
+ {
+ throw Exception("Unimplemented Configuration::data_coordinate_to_index");
+ /*
+ if( (!_receiver_on_emitter) )
+ {
+ throw Exception("Data space when an emitter is not also a receiver is not implemented yet");
+ }
+ std::array<unsigned int,2> index;
+
+ double DAngle = 2. * M_PI /((double) _np);
+ double Ds_e = DAngle * _emitter_skip * _rer;
+ double Ds_r = DAngle * _receiver_skip * _rer;
+ double Le0 = DAngle * _emitter_offset * _rer;
+ double Lr0 = DAngle * _receiver_offset * _rer;
+ index[0] = std::round( (xs-Le0) / Ds_e );
+ index[1] = std::round( (xr-Lr0) / Ds_r );
+
+ return index;
+ */
+ }
+
+ std::array<double,2> Configuration::index_to_data_coordinate(unsigned int s, unsigned int r) const
+ {
+ throw Exception("Unimplemented Configuration::index_to_data_coordinate");
+
+ /*
+ if( (!_receiver_on_emitter) )
+ {
+ throw Exception("Data space when an emitter is not also a receiver is not implemented yet");
+ }
+
+ std::array<double,2> index;
+
+ double DAngle = 2. * M_PI /((double) _np);
+ double Ds_e = DAngle * _emitter_skip * _rer;
+ double Ds_r = DAngle * _receiver_skip * _rer;
+ double Le0 = DAngle * _emitter_offset * _rer;
+ double Lr0 = DAngle * _receiver_offset * _rer;
+ index[0] = Le0 + ((double)s) * Ds_e;
+ index[1] = Lr0 + ((double)r) * Ds_r;
+
+ return index;*/
+ }
+
+ double Configuration::data_area() const
+ {
+ throw Exception("Unimplemented Configuration::data_area");
+ /*
+ double DAngle = 2. * M_PI /((double) _np);
+ double Ds_e = DAngle * _emitter_skip * _rer;
+ double Ds_r = DAngle * _receiver_skip * _rer;
+ double Le = (_ne-1) * Ds_e;
+ double Lr = (_nr-1) * Ds_r;
+
+ return Le * Lr;
+ */
+ };
+
+ double Configuration::datapoint_area() const
+ {
+ throw Exception("Unimplemented Configuration::datapoint_area");
+ /*
+ double DAngle = 2. * M_PI /((double) _np);
+ double Ds_e = DAngle * _emitter_skip * _rer;
+ double Ds_r = DAngle * _receiver_skip * _rer;
+
+ return Ds_e * Ds_r;*/
+ };
+
+ bool Configuration::data_coordinate_isValid(double xs,double xr) const
+ {
+ throw Exception("Unimplemented Configuration::data_coordinate_isValid");
+
+ /*
+ double DAngle = 2. * M_PI /((double) _np);
+ double Ds_e = DAngle * _emitter_skip * _rer;
+ double Ds_r = DAngle * _receiver_skip * _rer;
+ double Le0 = DAngle * _emitter_offset * _rer;
+ double Lr0 = DAngle * _receiver_offset * _rer;
+ double Le = (_ne-1) * Ds_e;
+ double Lr = (_nr-1) * Ds_r;
+
+ double eps = 1e-8;
+ gmsh::option::getNumber("Geometry.MatchMeshTolerance",eps);
+ if( (Le0-eps <= xs) && (xs <= Le0+Le+eps) && (Lr0-eps <= xr) && (xr <= Lr0+Lr+eps) )
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+
+ */
+ }
+
+ template<Physic T_Physic>
+ ModelMonoFunction green0_preconditioner(const Data<T_Physic>& dd, const WaveMultiField<T_Physic>& g, const Configuration* const config)
+ {
+ return autocorrelate<T_Physic>(g,config->emitter_idx())*autocorrelate<T_Physic>(g,config->receiver_idx());
+ }
+
+ template ModelMonoFunction green0_preconditioner<Physic::acoustic>(const Data<Physic::acoustic>&, const WaveMultiField<Physic::acoustic>&, const Configuration* const);
+ template ModelMonoFunction green0_preconditioner<Physic::electromagnetic>(const Data<Physic::electromagnetic>&, const WaveMultiField<Physic::electromagnetic>&, const Configuration* const);
+ template ModelMonoFunction green0_preconditioner<Physic::elastodynamic>(const Data<Physic::elastodynamic>&, const WaveMultiField<Physic::elastodynamic>&, const Configuration* const);
+
+} // namespace circular_acquisition
diff --git a/specific/configuration/flexible_acquisition.h b/specific/configuration/flexible_acquisition.h
new file mode 100644
index 0000000000000000000000000000000000000000..608b9e903ccfcd602714a9e29a12b0e78e874de7
--- /dev/null
+++ b/specific/configuration/flexible_acquisition.h
@@ -0,0 +1,121 @@
+#ifndef H_CONFIGURATION_FLEXIBLE_ACQUISITION
+#define H_CONFIGURATION_FLEXIBLE_ACQUISITION
+
+//GmshFEM Library
+#include "GmshFem.h"
+//#include "Function.h"
+
+//GmshFWI Library
+#include "../../common/configuration.h"
+#include "../../common/wave/element.h"
+#include "inclusion/inclusion.h"
+
+//Forward declaration
+template<Physic T_Physic>
+class Data;
+
+namespace flexible_acquisition
+{
+
+ class Configuration final : public ConfigurationInterface
+ {
+ private:
+
+ /**
+ * Emetters and receivers are stored together, and for each shot, one point is set as emitter,
+ * and some others (potentially including the emitter) are enabled as receivers
+ */
+ std::vector<std::pair<double, double>> _er_positions;
+
+ struct ShotConfig {
+ unsigned int emetterIdx;
+ std::vector<unsigned int> enabledReceivers;
+ };
+
+ std::vector<ShotConfig> _shotsConfigs;
+
+ /* Domain is a rectangle. "x" goes from 0 to L and "y" from -H to 0. Known part is y between 0 and -Hsub*/
+ double _H;
+ double _L;
+ double _Hsub;
+
+ // Mesh resolution
+ double _h;
+
+ // TODO
+ std::vector<unsigned int> _emitter_idx;
+ std::vector<unsigned int> _receiver_idx;
+
+ bool _areFilled;
+
+ // Subdomains: all in the parent class
+ /*
+ double _rer;
+ double _rext;
+
+ unsigned int _ne;
+ unsigned int _nr;
+ unsigned int _emitter_offset;
+ unsigned int _receiver_offset;
+ unsigned int _emitter_skip;
+ unsigned int _receiver_skip;
+ bool _receiver_on_emitter;
+
+ unsigned int _ni;
+
+ std::vector<const InclusionInterface*> _inclusion_geo;
+
+
+ UnknownRegion _unknown_region;
+ std::array<gmshfem::domain::Domain,2> _background;
+ std::array<gmshfem::domain::Domain,2> _inclusions;
+ std::vector<std::array<gmshfem::domain::Domain,2>> _inclusion;
+
+
+
+ std::vector<std::complex<double>> _mb;
+ std::vector<std::vector<std::complex<double>>> _mi;
+*/
+
+ std::array<gmshfem::domain::Domain, 2> _subsurface;
+ std::array<gmshfem::domain::Domain, 2> _supersurface;
+
+ std::vector<std::complex<double>> _m_super;
+ std::vector<std::complex<double>> _m_sub;
+
+ virtual void wave_mesh() const;
+ virtual void data_mesh() const;
+
+
+ void setupInitialModel(const ParametrizationInterface* const parametrization, const gmshfem::common::GmshFem& gmshFem);
+ public:
+ Configuration(std::string name, const ParametrizationInterface* const parametrization, const gmshfem::common::GmshFem& gmshFem);
+
+ const std::vector<unsigned int>& emitter_idx() const {return _emitter_idx;};
+ const std::vector<unsigned int>& receiver_idx() const {return _receiver_idx;};
+
+ bool areFilled() const {return _areFilled;};
+ unsigned int ner() const {return _np;};
+
+// virtual void mesh() const;
+
+ virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const;
+ virtual std::array<double,2> index_to_data_coordinate(unsigned int s, unsigned int r) const;
+ virtual bool data_coordinate_isValid(double xs,double xr) const;
+
+ virtual double area() const {return _H * _L;};
+ virtual double data_area() const override;
+ virtual double datapoint_area() const;
+ virtual double array() const {return 2. * (_L + _H);};
+ virtual double depth() const {return _H;};
+
+ virtual std::string wave_gmodel() const {return _name;};
+ virtual std::string model_gmodel() const {return _name;};
+ virtual std::string data_gmodel() const {return _name;};
+ };
+
+ template<Physic T_Physic>
+ ModelMonoFunction green0_preconditioner(const Data<T_Physic>& dd, const WaveMultiField<T_Physic>& g, const Configuration* const config);
+} // namespace soil
+
+#endif // H_CONFIGURATION_FLEXIBLE_ACQUISITION
diff --git a/specific/configuration/layeredvolume_acquisition.cpp b/specific/configuration/layeredvolume_acquisition.cpp
index a12e5a2541ac8e926c14ea25a28519d528944c29..a8399ea86b1e6050e08ee1623f49ee56fa64979d 100644
--- a/specific/configuration/layeredvolume_acquisition.cpp
+++ b/specific/configuration/layeredvolume_acquisition.cpp
@@ -480,26 +480,6 @@ namespace layeredvolume_acquisition
}
}
- void Configuration::mesh() const
- {
- msg::print << "Generate meshes" << msg::endl;
- gmsh::option::setNumber("General.Terminal", 1);
-
- gmodel::add(_name);
- if(!_remesh)
- {
- gmsh::open(_name+".msh");
- return;
- }
-
- wave_mesh();
- data_mesh();
-
- factory::synchronize();
- gmodel::mesh::generate();
-
- gmsh::write(_name+".msh");
- }
void Configuration::wave_mesh() const
{
std::vector<int> pl1(_layercoord1.size());
diff --git a/specific/configuration/layeredvolume_acquisition.h b/specific/configuration/layeredvolume_acquisition.h
index bb83d5d84cd3d6e585137eebe82b11686d24a47e..6618487b3f0b985f4c3579d8c9fe568e2c4b10d1 100644
--- a/specific/configuration/layeredvolume_acquisition.h
+++ b/specific/configuration/layeredvolume_acquisition.h
@@ -107,7 +107,6 @@ namespace layeredvolume_acquisition
unsigned int nye() const {return _nye;};
unsigned int nyr() const {return _nyr;};
- virtual void mesh() const;
virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const;
virtual std::array<double,2> index_to_data_coordinate(unsigned int s, unsigned int r) const;
diff --git a/specific/configuration/newConfiguration.h b/specific/configuration/newConfiguration.h
index 868b390d861218f40258c2360f816eb4fba6f420..9dd88c8292488df869a7405a5d8e987db91b1627 100644
--- a/specific/configuration/newConfiguration.h
+++ b/specific/configuration/newConfiguration.h
@@ -6,6 +6,7 @@
//Standard Library
#include <string>
+#include <memory>
//GmshFem Library
#include "GmshFem.h"
@@ -18,20 +19,22 @@
#include "circular_acquisition.h"
#include "line_acquisition.h"
#include "rectangular_acquisition.h"
+#include "flexible_acquisition.h"
-ConfigurationInterface* newConfiguration(std::string name, const ParametrizationInterface* const parametrization, const gmshfem::common::GmshFem& gmshFem)
+std::unique_ptr<ConfigurationInterface> newConfiguration(std::string name, const ParametrizationInterface* const parametrization, const gmshfem::common::GmshFem& gmshFem)
{
std::string configuration;
if(!gmshFem.userDefinedParameter(configuration, "configuration"))
{
throw gmshfem::common::Exception("Configuration type could not be found.");
}
- if(configuration=="surface_acquisition"){return new surface_acquisition::Configuration(name, parametrization, gmshFem);}
- else if(configuration=="volume_acquisition"){return new volume_acquisition::Configuration(name, parametrization, gmshFem);}
- else if(configuration=="layeredvolume_acquisition"){return new layeredvolume_acquisition::Configuration(name, parametrization, gmshFem);}
- else if(configuration=="circular_acquisition"){return new circular_acquisition::Configuration(name, parametrization, gmshFem);}
- else if(configuration=="line_acquisition"){return new line_acquisition::Configuration(name, parametrization, gmshFem);}
- else if(configuration=="rectangular_acquisition"){return new rectangular_acquisition::Configuration(name, parametrization, gmshFem);}
+ if(configuration=="surface_acquisition"){return std::make_unique< surface_acquisition::Configuration>(name, parametrization, gmshFem);}
+ else if(configuration=="volume_acquisition"){return std::make_unique< volume_acquisition::Configuration>(name, parametrization, gmshFem);}
+ else if(configuration=="layeredvolume_acquisition"){return std::make_unique< layeredvolume_acquisition::Configuration>(name, parametrization, gmshFem);}
+ else if(configuration=="circular_acquisition"){return std::make_unique< circular_acquisition::Configuration>(name, parametrization, gmshFem);}
+ else if(configuration=="line_acquisition"){return std::make_unique< line_acquisition::Configuration>(name, parametrization, gmshFem);}
+ else if(configuration=="rectangular_acquisition"){return std::make_unique< rectangular_acquisition::Configuration>(name, parametrization, gmshFem);}
+ else if(configuration=="flexible_acquisition"){return std::make_unique< flexible_acquisition::Configuration>(name, parametrization, gmshFem);}
else
{
throw gmshfem::common::Exception("Configuration" + configuration + " is not valid.");
diff --git a/specific/configuration/surface_acquisition.cpp b/specific/configuration/surface_acquisition.cpp
index 5f2a488ec3bb88a5fffd453fd370b8ba5e501f73..dac31c5fab66bf6f5e375d55a1d76ce5b8ea97d7 100644
--- a/specific/configuration/surface_acquisition.cpp
+++ b/specific/configuration/surface_acquisition.cpp
@@ -116,6 +116,8 @@ namespace surface_acquisition
_supersurface[Support::BLK] = Domain(2, 22);
_supersurface[Support::BND] = Domain(1, 12);
+ _named_domains["top_bnd"] = Domain("top_bnd");
+
_wave_omega[Support::BLK] = _subsurface[Support::BLK] | _supersurface[Support::BLK];
_wave_omega[Support::BND] = _subsurface[Support::BND] | _supersurface[Support::BND];
@@ -268,24 +270,7 @@ namespace surface_acquisition
}
}
- void Configuration::mesh() const
- {
- msg::print << "Generate meshes" << msg::endl;
- gmsh::option::setNumber("General.Terminal", 1);
- gmodel::add(_name);
- if(!_remesh)
- {
- gmsh::open(_name+".msh");
- return;
- }
- wave_mesh();
- data_mesh();
-
- factory::synchronize();
- gmodel::mesh::generate();
-
- gmsh::write(_name+".msh");
- }
+
void Configuration::wave_mesh() const
{
double Le0 = (_L - _Ler)/2.;
@@ -361,6 +346,7 @@ namespace surface_acquisition
lymin_reverse[n] = -lymin[n];
}
+ // Nomenclature: line from one y to another, left or right. All lines go up
int lHyminl = factory::addLine(pHl,pyminl);
int lymin0l = factory::addLine(pyminl,p0l);
int l0ymaxl = factory::addLine(p0l,pymaxl);
@@ -450,6 +436,10 @@ namespace surface_acquisition
gmodel::addPhysicalGroup(1, lsuper, 12);
gmodel::setPhysicalName(1, 12, "supersurface_bnd");
+ // Add another physical groupe for the top boundary, to represent e.g. reflection at the air-water interface
+ gmodel::addPhysicalGroup(1, lymax, 13);
+ gmodel::setPhysicalName(1, 13, "top_bnd");
+
gmodel::addPhysicalGroup(2, {s_sub}, 21);
gmodel::setPhysicalName(2, 21, "subsurface_vol");
gmodel::addPhysicalGroup(2, {s_supermin,s_supermax}, 22);
@@ -527,6 +517,10 @@ namespace surface_acquisition
gmodel::addPhysicalGroup(1, lsuper, 12);
gmodel::setPhysicalName(1, 12, "supersurface_bnd");
+ // Add another physical groupe for the top boundary, to represent e.g. reflection at the air-water interface
+ gmodel::addPhysicalGroup(1, lymax, 13);
+ gmodel::setPhysicalName(1, 13, "top_bnd");
+
gmodel::addPhysicalGroup(2, {s_sub}, 21);
gmodel::setPhysicalName(2, 21, "subsurface_vol");
gmodel::addPhysicalGroup(2, {s_supermax}, 22);
diff --git a/specific/configuration/surface_acquisition.h b/specific/configuration/surface_acquisition.h
index b4402359ac9663fd388b20b34a54a1db3a957062..00ad7d52aeee7904f124dfedcbaafd95c5fd50d1 100644
--- a/specific/configuration/surface_acquisition.h
+++ b/specific/configuration/surface_acquisition.h
@@ -59,7 +59,6 @@ namespace surface_acquisition
virtual void data_mesh() const;
public:
Configuration(std::string name, const ParametrizationInterface* const parametrization, const gmshfem::common::GmshFem& gmshFem);
- virtual void mesh() const;
unsigned int ne() const {return _ne;};
unsigned int nr() const {return _nr;};
diff --git a/specific/configuration/volume_acquisition.cpp b/specific/configuration/volume_acquisition.cpp
index a8fe8673ab9d844c87346d7f613fc6c041b9987e..082a57ef1910b98c32c1797e79461ec8d5dbb5f3 100644
--- a/specific/configuration/volume_acquisition.cpp
+++ b/specific/configuration/volume_acquisition.cpp
@@ -344,26 +344,6 @@ namespace volume_acquisition
}
}
- void Configuration::mesh() const
- {
- msg::print << "Generate meshes" << msg::endl;
- gmsh::option::setNumber("General.Terminal", 1);
-
- gmodel::add(_name);
- if(!_remesh)
- {
- gmsh::open(_name+".msh");
- return;
- }
-
- wave_mesh();
- data_mesh();
-
- factory::synchronize();
- gmodel::mesh::generate();
-
- gmsh::write(_name+".msh");
- }
void Configuration::wave_mesh() const
{
int pb1 = factory::addPoint(0., 0., 0., _h);
diff --git a/specific/configuration/volume_acquisition.h b/specific/configuration/volume_acquisition.h
index faf51aa50f6c81a2b8421e546d1827ee9a26640d..302cc85e7e33d30a53a3334e673b42726cfbb2ee 100644
--- a/specific/configuration/volume_acquisition.h
+++ b/specific/configuration/volume_acquisition.h
@@ -92,7 +92,6 @@ namespace volume_acquisition
unsigned int nye() const {return _nye;};
unsigned int nyr() const {return _nyr;};
- virtual void mesh() const;
virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const;
virtual std::array<double,2> index_to_data_coordinate(unsigned int s, unsigned int r) const;
diff --git a/specific/configuration/yaml_interface.cpp b/specific/configuration/yaml_interface.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a7baa4dea954bb6cd5a1444367303c06cae51766
--- /dev/null
+++ b/specific/configuration/yaml_interface.cpp
@@ -0,0 +1,34 @@
+#include "GmshFem.h"
+#include "yaml_interface.h"
+
+ShotsConfigurationYAML::ShotsConfigurationYAML(const std::string &path)
+{
+#ifndef ENABLE_YAML_CPP
+ throw gmshfem::common::Exception("Cannot use Flexible Acquisition without YAML-CPP.");
+#else
+ YAML::Node config = YAML::LoadFile(path);
+
+ for (const auto &coord : config["coordinates"])
+ {
+ double x = coord[0].as<double>();
+ double y = coord[1].as<double>();
+ _points.emplace_back(x, y);
+ }
+
+ for (const auto &shot : config["shots"])
+ {
+ std::vector<unsigned> emitters, receivers;
+ for (const auto &emitter : shot["emitters"])
+ {
+ emitters.push_back(emitter.as<int>());
+ }
+ for (const auto &receiver : shot["receivers"])
+ {
+ receivers.push_back(receiver.as<int>());
+ }
+
+ _emitters.push_back(std::move(emitters));
+ _receivers.push_back(std::move(receivers));
+ }
+#endif
+}
\ No newline at end of file
diff --git a/specific/configuration/yaml_interface.h b/specific/configuration/yaml_interface.h
new file mode 100644
index 0000000000000000000000000000000000000000..23afb3fcdbde94a5e01e4d678bbba79a1e72f0d8
--- /dev/null
+++ b/specific/configuration/yaml_interface.h
@@ -0,0 +1,77 @@
+#ifndef H_YAML_INTERFACE_SHOTS
+#define H_YAML_INTERFACE_SHOTS
+
+#include <vector>
+#include <string>
+#include <utility>
+#include "yaml-cpp/yaml.h"
+// TODO: refactor (abstract + impl)
+
+class ShotsConfigurationInterface //(abstract)
+{
+ public:
+ using Coordinates = std::pair<double, double>;
+
+ /**
+ * Describe the number of shots
+ */
+ virtual unsigned numShots() const = 0;
+
+ /**
+ * Describe the number of points
+ */
+ virtual unsigned numPoints() const = 0;
+
+ /**
+ * emitters(s) is the list of enabled emitters in shot "s".
+ * Array must be of size numShots()
+ */
+ virtual const std::vector<unsigned> &emitters(unsigned shot) const = 0;
+
+ /**
+ * receivers(s) is the list of enabled receivers in shot "s".
+ * Array must be of size numShots()
+ */
+ virtual const std::vector<unsigned> &receivers(unsigned shot) const = 0;
+
+ /**
+ * List of all points
+ * Array must be of size numPoints()
+ */
+ virtual const std::vector<Coordinates> &points() const = 0;
+};
+
+class ShotsConfigurationYAML : public ShotsConfigurationInterface {
+ private:
+ std::vector<std::vector<unsigned>> _emitters, _receivers;
+ std::vector<ShotsConfigurationInterface::Coordinates> _points;
+
+ public:
+
+ ShotsConfigurationYAML(const std::string& path);
+
+ virtual unsigned numShots() const override {
+ return _emitters.size();
+ }
+
+ virtual unsigned numPoints() const override {
+ return _points.size();
+ }
+
+ virtual const std::vector<unsigned> &emitters(unsigned shot) const override {
+ return _emitters.at(shot);
+ }
+
+ virtual const std::vector<unsigned> &receivers(unsigned shot) const override {
+ return _receivers.at(shot);
+ }
+
+ virtual const std::vector<Coordinates> &points() const override {
+ return _points;
+ }
+
+
+
+};
+
+#endif
\ No newline at end of file
diff --git a/specific/data/element.cpp b/specific/data/element.cpp
index 26a948ad32be0484bae5680035e66487a479bda4..6bef6003784383cd330db1fec1bd21b59c0fbedf 100644
--- a/specific/data/element.cpp
+++ b/specific/data/element.cpp
@@ -12,6 +12,13 @@
using namespace gmshfem;
using namespace gmshfem::common;
+template<>
+typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::acoustic>::degree >::Object s_ones<Physic::acoustic>(){return 1.;}
+template<>
+typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::electromagnetic>::degree >::Object s_ones<Physic::electromagnetic>(){return PointData<Physic::electromagnetic>::Object(1.,1.,1.);}
+template<>
+typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::elastodynamic>::degree >::Object s_ones<Physic::elastodynamic>(){return PointData<Physic::elastodynamic>::Object(1.,1.,1.);}
+
/* acoustic */
template<>
double norm2<Physic::acoustic>(const typename PointData<Physic::acoustic>::Object& d){return std::norm(d);}
@@ -27,8 +34,6 @@ template<>
typename PointData<Physic::acoustic>::Object wise_divide<Physic::acoustic>(const typename PointData<Physic::acoustic>::Object& d1,const typename PointData<Physic::acoustic>::Object& d2){return d1/d2;}
template<>
typename PointData<Physic::acoustic>::Object wise_multiply<Physic::acoustic>(const typename PointData<Physic::acoustic>::Object& d1,const typename PointData<Physic::acoustic>::Object& d2){return d1*d2;}
-template<>
-typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::acoustic>::degree >::Object s_ones<Physic::acoustic>(){return 1.;}
/* electromagnetic */
template<>
@@ -48,8 +53,6 @@ template<>
typename PointData<Physic::electromagnetic>::Object wise_divide<Physic::electromagnetic>(const typename PointData<Physic::electromagnetic>::Object& d1,const typename PointData<Physic::electromagnetic>::Object& d2){return d1.array() / d2.array() ;}
template<>
typename PointData<Physic::electromagnetic>::Object wise_multiply<Physic::electromagnetic>(const typename PointData<Physic::electromagnetic>::Object& d1,const typename PointData<Physic::electromagnetic>::Object& d2){return d1.array() * d2.array() ;}
-template<>
-typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::electromagnetic>::degree >::Object s_ones<Physic::electromagnetic>(){return PointData<Physic::electromagnetic>::Object(1.,1.,1.);}
/* elastodynamic */
template<>
@@ -69,5 +72,3 @@ template<>
typename PointData<Physic::elastodynamic>::Object wise_divide<Physic::elastodynamic>(const typename PointData<Physic::elastodynamic>::Object& d1,const typename PointData<Physic::elastodynamic>::Object& d2){return d1.array() / d2.array() ;}
template<>
typename PointData<Physic::elastodynamic>::Object wise_multiply<Physic::elastodynamic>(const typename PointData<Physic::elastodynamic>::Object& d1,const typename PointData<Physic::elastodynamic>::Object& d2){return d1.array() * d2.array() ;}
-template<>
-typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::elastodynamic>::degree >::Object s_ones<Physic::elastodynamic>(){return PointData<Physic::elastodynamic>::Object(1.,1.,1.);}
diff --git a/specific/data/objective/newObjective.h b/specific/data/objective/newObjective.h
index 67c43e8574ef39c53f2c51de9b7987b82c8639e9..193f9966588f32ffc454340dfcb281d9cd61d94e 100644
--- a/specific/data/objective/newObjective.h
+++ b/specific/data/objective/newObjective.h
@@ -6,6 +6,7 @@
//Standard Library
#include <string>
+#include <memory>
//GmshFEM Library
#include "Exception.h"
@@ -15,16 +16,16 @@
#include "phase.h"
template<Physic T_Physic>
-ObjectiveInterface<T_Physic>* newObjective(const Data<T_Physic>& d0, const ConfigurationInterface* const config, const gmshfem::common::GmshFem& gmshFem,std::string suffix = "")
+std::unique_ptr<ObjectiveInterface<T_Physic>> newObjective(const Data<T_Physic>& d0, const ConfigurationInterface* const config, const gmshfem::common::GmshFem& gmshFem,std::string suffix = "")
{
std::string objective;
if(!gmshFem.userDefinedParameter(objective, "objective"))
{
throw gmshfem::common::Exception("Objective type could not be found.");
}
- if(objective=="l2distance"){return new l2distance::Objective<T_Physic>(d0,gmshFem);}
- else if(objective=="conventional_phase"){return new conventional_phase::Objective<T_Physic>(d0,gmshFem);}
- else if(objective=="logarithmic_phase"){return new logarithmic_phase::Objective<T_Physic>(d0,gmshFem);}
+ if(objective=="l2distance"){return std::make_unique< l2distance::Objective<T_Physic> >(d0,gmshFem);}
+ else if(objective=="conventional_phase"){return std::make_unique< conventional_phase::Objective<T_Physic> >(d0,gmshFem);}
+ else if(objective=="logarithmic_phase"){return std::make_unique< logarithmic_phase::Objective<T_Physic> >(d0,gmshFem);}
else
{
throw gmshfem::common::Exception("Objective " + objective + " is not valid.");
diff --git a/specific/model/innerproduct/newInnerProduct.cpp b/specific/model/innerproduct/newInnerProduct.cpp
index b7b8f38d47039af1bf65f54d69dc99d9f1a63d7d..651af62ad665b8cadea536c17d212599507cb4ed 100644
--- a/specific/model/innerproduct/newInnerProduct.cpp
+++ b/specific/model/innerproduct/newInnerProduct.cpp
@@ -11,15 +11,15 @@
#include "sobolev.h"
template<class M>
-InnerProductInterface* newInnerProduct(const ConfigurationInterface* const config, const M& m_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix)
+std::unique_ptr<InnerProductInterface> newInnerProduct(const ConfigurationInterface* const config, const M& m_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix)
{
std::string innerproduct;
if(!gmshFem.userDefinedParameter(innerproduct, "innerproduct"))
{
throw gmshfem::common::Exception("Inner product type could not be found.");
}
- if(innerproduct=="sobolev_l2"){return new sobolev_l2::InnerProduct(config,m_discret,gmshFem, suffix);}
- else if(innerproduct=="sobolev_df"){return new sobolev_df::InnerProduct(config,m_discret,gmshFem, suffix);}
+ if(innerproduct=="sobolev_l2"){return std::make_unique< sobolev_l2::InnerProduct>(config,m_discret,gmshFem, suffix);}
+ else if(innerproduct=="sobolev_df"){return std::make_unique< sobolev_df::InnerProduct>(config,m_discret,gmshFem, suffix);}
else
{
throw gmshfem::common::Exception("Inner product " + innerproduct + " is not valid.");
@@ -27,5 +27,5 @@ InnerProductInterface* newInnerProduct(const ConfigurationInterface* const confi
}
}
-template InnerProductInterface* newInnerProduct<model::Discretization>(const ConfigurationInterface* const, const model::Discretization&, const gmshfem::common::GmshFem&, std::string);
-template InnerProductInterface* newInnerProduct<ModelField>(const ConfigurationInterface* const, const ModelField&, const gmshfem::common::GmshFem&, std::string);
+template std::unique_ptr<InnerProductInterface> newInnerProduct<model::Discretization>(const ConfigurationInterface* const, const model::Discretization&, const gmshfem::common::GmshFem&, std::string);
+template std::unique_ptr<InnerProductInterface> newInnerProduct<ModelField>(const ConfigurationInterface* const, const ModelField&, const gmshfem::common::GmshFem&, std::string);
diff --git a/specific/model/innerproduct/newInnerProduct.h b/specific/model/innerproduct/newInnerProduct.h
index c968ffac047540e6b04db4d4d3f879362d43c530..dc77005f3af143108dccfa86a77cdb6523e7f066 100644
--- a/specific/model/innerproduct/newInnerProduct.h
+++ b/specific/model/innerproduct/newInnerProduct.h
@@ -5,12 +5,13 @@
#define H_SPECIFIC_MODEL_NEWINNERPRODUCT
//Standard Library
+#include <memory>
//GmshFEM Library
#include "GmshFem.h"
//GmshFWI Library
#include "../../../common/model/innerproduct/innerproduct.h"
template<class M>
-InnerProductInterface* newInnerProduct(const ConfigurationInterface* const config, const M& m_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix);
+std::unique_ptr<InnerProductInterface> newInnerProduct(const ConfigurationInterface* const config, const M& m_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix);
#endif //H_SPECIFIC_MODEL_NEWINNERPRODUCT
diff --git a/specific/optimization/globalminimumsearch/sequential.cpp b/specific/optimization/globalminimumsearch/sequential.cpp
index 500aa5fd5d05e878ad8847ffd408302f27c404f9..1aa559a93f94b712b680c2664fca96b6eb028cea 100644
--- a/specific/optimization/globalminimumsearch/sequential.cpp
+++ b/specific/optimization/globalminimumsearch/sequential.cpp
@@ -15,8 +15,13 @@
#include "../../optimization/linesearch/newLineSearch.h"
#include "../../optimization/localminimumsearch/newLocalMinimumSearch.h"
+// STL
+#include <memory>
+
using namespace gmshfem;
using namespace gmshfem::common;
+using std::unique_ptr;
+using std::make_unique;
namespace sequential
{
@@ -143,12 +148,11 @@ namespace sequential
throw Exception("Data and equations do not have the same size.");
}
- ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(d0,_config,_gmshFem,_suffix);
- InnerProductInterface* const innerproduct = newInnerProduct(_config,*m,_gmshFem,_suffix);
+ auto objective = newObjective<T_Physic>(d0,_config,_gmshFem,_suffix);
+ auto const innerproduct = newInnerProduct(_config,*m,_gmshFem,_suffix);
RegularizationInterface* const regularization = newRegularization(_config,*m,_gmshFem,_suffix);
- StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(_config,innerproduct,regularization,d0.frequency(),equation,objective);
- FunctionalInterface* const functional = statefunctional;
- functional->setModel(*m);
+ auto statefunctional = make_unique<StateFunctional<T_Physic>>(_config,innerproduct.get(),regularization,d0.frequency(),equation,objective.get());
+ statefunctional->setModel(*m);
LocalMinimumSearchInterface* const minimumsearch = newLocalMinimumSearch(_gmshFem);
DescentSearchInterface* const descentsearch = newDescentSearch(_gmshFem);
@@ -157,7 +161,7 @@ namespace sequential
std::vector<GlobalMinimumSearchHistoryLine2> historyline;
double js_1 =0.;
- double js = functional->performance();
+ double js = statefunctional->performance();
historyline.emplace_back(js,0.,0.,statefunctional->performance_objective(),statefunctional->performance_regularization());
double rel_djs = 0.;
bool success = false;
@@ -184,8 +188,8 @@ namespace sequential
}
minimumsearch->name(_name+"g"+_suffix+"s"+suffix1);
js_1=js;
- (*minimumsearch)(m, functional, descentsearch, linesearch);
- js = functional->performance();
+ (*minimumsearch)(m, statefunctional.get(), descentsearch, linesearch);
+ js = statefunctional->performance();
m->write(_name+"_mg"+_suffix+"s"+suffix1);
//statefunctional->write_data(Type::FS,_name+"_dg"+_suffix+"s"+suffix1,"pos");
rel_djs = 2. * (js_1-js) / (js_1+js);
@@ -207,10 +211,7 @@ namespace sequential
delete descentsearch;
delete linesearch;
- delete functional;
- delete innerproduct;
delete regularization;
- delete objective;
msg::unindent();
}
diff --git a/specific/optimization/localminimumsearch/classic.cpp b/specific/optimization/localminimumsearch/classic.cpp
index b0498aa39cffde89afbb47458457a1cf4ad8ab06..47f6e4a7838d001c181cd95ca629a22c4771dff2 100644
--- a/specific/optimization/localminimumsearch/classic.cpp
+++ b/specific/optimization/localminimumsearch/classic.cpp
@@ -82,6 +82,13 @@ namespace classic
msg::print << "Search minimum (classic)" << msg::endl;
msg::indent();
+ if (!linesearch) {
+ throw Exception("No linesearch provided.");
+ }
+ if (!descentsearch) {
+ throw Exception("No descentsearch provided.");
+ }
+
descentsearch->initialize(functional);
std::vector<LocalMinimumSearchHistoryLine2> historyline;
diff --git a/specific/wave/discretization.cpp b/specific/wave/discretization.cpp
index 3282e8ce822dec5e0c5d6b446aad9641fbbdc329..f27798f786830f5e824989aecc0c9a6b25755810 100644
--- a/specific/wave/discretization.cpp
+++ b/specific/wave/discretization.cpp
@@ -15,37 +15,6 @@ using namespace gmshfem::field;
namespace wave
{
- /*
- * class Discretization
- */
- template<Physic T_physic>
- Discretization<T_physic>::Discretization(const GmshFem& gmshFem, std::string suffix)
- {
- msg::print << "Read wave discretization parameters" << msg::endl;
- std::string str_buffer;
- if(!(
- (
- gmshFem.userDefinedParameter(_functionSpaceDegree, "wave_FunctionSpaceDegree")
- ||
- gmshFem.userDefinedParameter(_functionSpaceDegree, "wave_FunctionSpaceDegree"+suffix)
- )
- &&
- (
- gmshFem.userDefinedParameter(str_buffer, "wave_FunctionSpaceType")
- ||
- gmshFem.userDefinedParameter(str_buffer, "wave_FunctionSpaceType"+suffix)
- )
- ))
- {
- throw common::Exception("A discretization parameter could not be found.");
- }
- _functionSpaceType = to_functionSpaceType<T_physic>(str_buffer);
- }
-
- template class Discretization<Physic::acoustic>;
- template class Discretization<Physic::electromagnetic>;
- template class Discretization<Physic::elastodynamic>;
-
/*
* to_functionSpaceType
*/
@@ -81,4 +50,36 @@ namespace wave
else { throw common::Exception("Invalid function space type for elastodynamics."); }
return out;
};
+
+ /*
+ * class Discretization
+ */
+ template<Physic T_physic>
+ Discretization<T_physic>::Discretization(const GmshFem& gmshFem, std::string suffix)
+ {
+ msg::print << "Read wave discretization parameters" << msg::endl;
+ std::string str_buffer;
+ if(!(
+ (
+ gmshFem.userDefinedParameter(_functionSpaceDegree, "wave_FunctionSpaceDegree")
+ ||
+ gmshFem.userDefinedParameter(_functionSpaceDegree, "wave_FunctionSpaceDegree"+suffix)
+ )
+ &&
+ (
+ gmshFem.userDefinedParameter(str_buffer, "wave_FunctionSpaceType")
+ ||
+ gmshFem.userDefinedParameter(str_buffer, "wave_FunctionSpaceType"+suffix)
+ )
+ ))
+ {
+ throw common::Exception("A discretization parameter could not be found.");
+ }
+ _functionSpaceType = to_functionSpaceType<T_physic>(str_buffer);
+ }
+
+ template class Discretization<Physic::acoustic>;
+ template class Discretization<Physic::electromagnetic>;
+ template class Discretization<Physic::elastodynamic>;
+
} // namespace wave
diff --git a/specific/wave/equation/simplewave.cpp b/specific/wave/equation/simplewave.cpp
index 56eba61d476dcbbc34c41ac14c936137e2e51900..a4ad2680d22cdf3e3e2b525030bedc7aa2ae1ea2 100644
--- a/specific/wave/equation/simplewave.cpp
+++ b/specific/wave/equation/simplewave.cpp
@@ -75,6 +75,21 @@ namespace simplewave
{
throw Exception("Impossible to compute boundary terms with real impedance: sensitivities are wrong.");
}
+
+ if (gmshFem.userDefinedParameter(_topBC, "top_bc")) {
+ if (_topBC != "Absorbing" && _topBC != "Dirichlet" && _topBC != "Neumann") {
+ throw Exception("Invalid value of parameter top_bc: " + _topBC);
+ }
+ }
+ else {
+ msg::warning << "SimpleWave without specified top_bc defaults to Absorbing." << msg::endl;
+ }
+
+ _useSecondOrderABC = false;
+ int intArg;
+ gmshFem.userDefinedParameter(intArg, "useSecondOrderABC");
+ _useSecondOrderABC = intArg != 0;
+ msg::warning << "Using Second Order ABC : " << _useSecondOrderABC << msg::endl;
}
bool Equation::modelIsObsolete()
@@ -107,16 +122,40 @@ namespace simplewave
}
void Equation::setLHS(const ModelStateEvaluator& ms)
{
+ auto top = _config->named_domains().find("top_bnd");
+
+
+ if (_topBC == "Dirichlet" || _topBC == "Neumann")
+ {
+ if (top == _config->named_domains().end())
+ throw Exception("Asked for specific BCs on top but this entity isn't defined.");
+ }
+
+ if (_topBC == "Dirichlet")
+ {
+ _v.addConstraint(top->second, 0.0);
+ }
+
_formulation.integral(-1.*grad(dof(_v)), grad(tf(_v)), _config->wave_omega(Support::BLK),integrationType(2*_w_degree));
_formulation.integral(_pulsation * _pulsation * m0 * dof(_v), tf(_v), _config->model_known(Support::BLK), integrationType(_integrationDegreeBlk));
+
+ auto absorbingPartKnown = _config->model_known(Support::BND);
+ auto absorbingPartUnknown = _config->model_unknown(Support::BND);
+ if (_topBC == "Neumann") {
+ absorbingPartKnown &= ~(top->second);
+ absorbingPartUnknown &= ~(top->second);
+ }
if(_realImpedance)
{
- _formulation.integral(-1. * im * _pulsation * sqrt(real(m0)) * dof(_v), tf(_v), _config->model_known(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ _formulation.integral(-1. * im * _pulsation * sqrt(real(m0)) * dof(_v), tf(_v), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ if (_useSecondOrderABC) _formulation.integral(1. * im / (2.0 * _pulsation * sqrt(real(m0))) * grad(dof(_v)), grad(tf(_v)), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+
}
else
{
- _formulation.integral(-1. * im * _pulsation * sqrt(m0) * dof(_v), tf(_v), _config->model_known(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ _formulation.integral(-1. * im * _pulsation * sqrt(m0) * dof(_v), tf(_v), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ if (_useSecondOrderABC) _formulation.integral(1. * im / (2.0 * _pulsation * sqrt(m0)) * grad(dof(_v)), grad(tf(_v)), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
if(!_config->model_unknown(Support::BLK).isEmpty())
{
@@ -126,11 +165,13 @@ namespace simplewave
{
if(_realImpedance)
{
- _formulation.integral(-1. * im * _pulsation * sqrt(real(m)) * dof(_v), tf(_v), _config->model_unknown(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ _formulation.integral(-1. * im * _pulsation * sqrt(real(m)) * dof(_v), tf(_v), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ if (_useSecondOrderABC) _formulation.integral(1. * im / (2. * _pulsation * sqrt(real(m))) * grad(dof(_v)), grad(tf(_v)), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
else
{
- _formulation.integral(-1. * im * _pulsation * sqrt(m) * dof(_v), tf(_v), _config->model_unknown(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ _formulation.integral(-1. * im * _pulsation * sqrt(m) * dof(_v), tf(_v), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ if (_useSecondOrderABC) _formulation.integral(1. * im / (2. * _pulsation * sqrt(m)) * grad(dof(_v)), grad(tf(_v)), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
}
}
diff --git a/specific/wave/equation/simplewave.h b/specific/wave/equation/simplewave.h
index 9e34a5815a36b976b719e5d0557ad7a401109621..2122ab80ac28f86bc43ba9ad74ee09194804bf08 100644
--- a/specific/wave/equation/simplewave.h
+++ b/specific/wave/equation/simplewave.h
@@ -24,6 +24,8 @@ namespace simplewave
bool _realImpedance;
double _const_preconditioner;
using EquationInterface<Physic::acoustic>::_boundary;
+ std::string _topBC;
+ bool _useSecondOrderABC;
public:
Equation(unsigned int f_idx,double pulsation,const ConfigurationInterface* const config, const wave::Discretization<Physic::acoustic>& w_discret,const gmshfem::common::GmshFem& gmshFem, std::string suffix = "");
diff --git a/synthetics.cpp b/synthetics.cpp
index 29bef54ecb1fa5b1a839e60bd9d534a26aab8a24..8f1a4a422512989adf8fa556bf8da2694f83901e 100644
--- a/synthetics.cpp
+++ b/synthetics.cpp
@@ -24,9 +24,9 @@ int synthetics(const GmshFem& gmshFem)
const ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
- ConfigurationInterface* configuration = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
- ModelState ms(configuration);
+ ModelState ms(configuration.get());
unsigned int n_freq = 1;
gmshFem.userDefinedParameter(n_freq, "n_freq");
@@ -57,10 +57,10 @@ int synthetics(const GmshFem& gmshFem)
msg::indent();
wave::Discretization<T_Physic> w_discret(gmshFem,suffix);
- EquationInterface<T_Physic>* const pequation = new ParametrizedEquation<T_Physic>(parametrization,0,2.*M_PI*frequency,configuration,w_discret,gmshFem,suffix);
+ EquationInterface<T_Physic>* const pequation = new ParametrizedEquation<T_Physic>(parametrization,0,2.*M_PI*frequency,configuration.get(),w_discret,gmshFem,suffix);
- WaveUpdater<T_Physic> wu(configuration,nullptr,pequation);
- DataUpdater<T_Physic> du(std::vector<double>(1,frequency),configuration,&wu,nullptr);
+ WaveUpdater<T_Physic> wu(configuration.get(),nullptr,pequation);
+ DataUpdater<T_Physic> du(std::vector<double>(1,frequency),configuration.get(),&wu,nullptr);
std::array<bool,5> dataNeedUpToDate = {true,false,false,false,false};
du.get(dataNeedUpToDate,ms).write(Type::FS,name+"_data"+suffix);
@@ -97,7 +97,6 @@ int synthetics(const GmshFem& gmshFem)
msg::unindent();
}
msg::unindent();
- delete configuration;
delete parametrization;
return 0;
diff --git a/tests/test.cpp b/tests/test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3a865294a01664b6f63317241c7bbd6c11c8c7f0
--- /dev/null
+++ b/tests/test.cpp
@@ -0,0 +1,109 @@
+#define BOOST_TEST_MODULE fwi_tests
+#include <boost/test/unit_test.hpp>
+
+#include <cstring>
+
+#include <gmsh.h>
+#include <gmshfem/Formulation.h>
+#include <gmshfem/GmshFem.h>
+#include "../common/wave/element.h"
+#include "../common/data/element.h"
+#include "../specific/data/element.h"
+
+#include "dummyAcquisition.h"
+
+#define INIT_GMSHFEM() \
+ namespace geo = gmsh::model::geo; \
+ namespace model = gmsh::model; \
+ namespace mesh = gmsh::model::mesh; \
+ using namespace gmshfem; \
+ using namespace gmshfem::common; \
+ using namespace gmshfem::problem; \
+ using namespace gmshfem::domain; \
+ using namespace gmshfem::field; \
+ using namespace gmshfem::function; \
+ using namespace gmshfem::post; \
+ using namespace gmshfem::equation; \
+ gmshfem::common::GmshFem fem(0, nullptr);
+
+BOOST_AUTO_TEST_SUITE(WaveMultiFieldTest);
+
+BOOST_AUTO_TEST_CASE(TestHomoegenuousDirichletWMF)
+{
+ INIT_GMSHFEM();
+ fem.getArgsManager()->addUserDefinedParameter("wave_FunctionSpaceDegree", "2");
+ fem.getArgsManager()->addUserDefinedParameter("wave_FunctionSpaceType", "HierarchicalH1");
+
+ double lc = 0.1;
+ model::add("Poisson");
+
+ geo::addPoint(0.0, 0.0, 0.0, lc, 1);
+ geo::addPoint(1.0, 0.0, 0.0, lc, 2);
+ geo::addPoint(1.0, 1.0, 0.0, lc, 3);
+ geo::addPoint(0.0, 1.0, 0.0, lc, 4);
+
+ geo::addLine(1, 2);
+ geo::addLine(2, 3);
+ geo::addLine(3, 4);
+ geo::addLine(4, 1);
+
+ geo::addCurveLoop({1,2,3,4});
+ geo::addPlaneSurface({1});
+ geo::synchronize();
+ mesh::generate();
+
+ model::addPhysicalGroup(1, {4}, 1);
+ model::setPhysicalName(1, 1, "gammaLeft");
+
+ model::addPhysicalGroup(2, {1}, 1);
+ model::setPhysicalName(2, 1, "omega");
+
+ Formulation<std::complex<double>> formulation("Poisson");
+ Domain omega("omega");
+ Domain gammaLeft("gammaLeft");
+
+ // Allocate an Auxiliary Field that will serve as proxy to write in the WMF
+ WaveAuxilaryField<Physic::acoustic> v("v", omega, "Poisson", wave::Discretization<Physic::acoustic>(fem));
+ v.addConstraint(gammaLeft, 0.); // Apply Dirichlet *before* preprocessing, to get correct DOF dictionary
+
+ // Assemble LHS to get the correct number of DOFs
+ formulation.integral(grad(dof(v)), grad(tf(v)), omega, "Gauss6");
+ formulation.pre();
+
+ WaveMultiField<Physic::acoustic> wmf(2, "v_multi");
+ wmf.initializeField(v);
+ v.setField(&wmf);
+
+ // First solve
+ v.setIndex(0);
+ formulation.integral(-1., tf(v), omega, "Gauss6");
+ formulation.assemble();
+ formulation.solve(true);
+
+ // Second solve
+ formulation.setRHSToZero();
+
+ v.setIndex(1);
+ formulation.removeTerms();
+ formulation.integral(-2., tf(v), omega, "Gauss6");
+ formulation.assemble();
+ formulation.solve(true);
+
+ // Check analytical solution to tolerance.
+ double f = 1;
+ auto expected = f * x<std::complex<double>>() * (1. - 0.5 * x<std::complex<double>>());
+
+ // Quadratic solution -> should be exact on a coarse mesh.
+ auto integralV = abs(integrate(v, omega, "Gauss4"));
+ auto integralErr1 = abs(integrate(expected - wmf[0], omega, "Gauss4"));
+ auto integralErr2 = abs(integrate(2. * expected - wmf[1], omega, "Gauss4"));
+ BOOST_TEST(integralV < 1e-8, "AuxiliaryWaveField should be empty. Integral of absolute value is " << integralV);
+ BOOST_TEST(integralErr1 < 1e-8, "First WMF has the wrong value. Error norm: " << integralErr1);
+ BOOST_TEST(integralErr2 < 1e-8, "Second WMF has the wrong value. Error norm: " << integralErr2);
+}
+BOOST_AUTO_TEST_CASE(dummy_test)
+{
+ BOOST_CHECK(1 + 1 == 2);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
diff --git a/time_synthetics_wave.cpp b/time_synthetics_wave.cpp
index 49c571c986e34416a039c45b6fd1fb82455f2c3f..590f427308b158927701c2fcbbed6ab3a00b514a 100644
--- a/time_synthetics_wave.cpp
+++ b/time_synthetics_wave.cpp
@@ -94,7 +94,7 @@ int main(int argc, char **argv)
const ParametrizationInterface* const parametrization = newParametrization<Physic::acoustic>(gmshFem);
- ConfigurationInterface* config = newConfiguration(name, parametrization, gmshFem);
+ auto configuration = newConfiguration(name, parametrization, gmshFem);
ModelState ms(config);