Permalink
Browse files

Merge branch 'dev' of https://github.com/gimli-org/gimli into dev

  • Loading branch information...
florian-wagner committed Aug 4, 2018
2 parents eb92741 + 8e199ef commit 11b53a8bd7973cd856a974d5a1180a20b22f74c3
View
@@ -18,6 +18,15 @@ endif()
project(libgimli)
# Check if conda package is created
if(DEFINED ENV{CONDA_BUILD})
message(STATUS "Conda package is being created.")
set(CONDA_BUILD TRUE)
set(Boost_INCLUDE_DIR "${CMAKE_PREFIX_PATH}/include")
else()
set(CONDA_BUILD FALSE)
endif()
if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
set(CMAKE_CXX_FLAGS_RELEASE "-O2 -pipe -ansi -Wall -Wno-long-long -Wno-unused-result -Wno-unused-variable")
@@ -26,37 +35,34 @@ if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" ST
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -p -O2 -ansi -pedantic -fno-omit-frame-pointer -Wall")
if (WIN32)
if (WIN32 OR CONDA_BUILD)
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -Wno-deprecated-declarations")
endif()
if(NOT WIN32 AND ASAN)
if (NOT WIN32 AND ASAN)
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fsanitize=address")
endif()
set(CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS "-Os")
# if (NOT WIN32)
# add_definitions(-std=c++0x)
add_definitions(-std=c++11)
#endif()
#add_definitions(-std=gnu++0x)
#-Wl,--no-undefined
#set(CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS "-Os -static-libgcc -Wl,-O2 -Wl,--as-needed -Wl,--sort-common")
#set(CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS "-Os -static-libgcc")
# set(CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS "-Os -Wl,--no-undefined -static-libgcc -Wl,-O2 -Wl,--as-needed -Wl,--sort-common")
# set(CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS "${CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS} -lbfd")
# add_definitions(-std=gnu++0x)
# -Wl,--no-undefined
# set(CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS "-Os -static-libgcc -Wl,-O2 -Wl,--as-needed -Wl,--sort-common")
# set(CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS "-Os -static-libgcc")
# set(CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS "-Os -Wl,--no-undefined -static-libgcc -Wl,-O2 -Wl,--as-needed -Wl,--sort-common")
# set(CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS "${CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS} -lbfd")
if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -Wno-overloaded-virtual")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -Qunused-arguments")
if (APPLE)
#If it's libc++ and you're on <= 10.8, you need to compile with clang++ -stdlib=libc++. If it's libstdc++ and you're on 10.9 or later, you need to compile with clang++ -stdlib=libstdc++.
#add_compile_options(-stdlib=libstdc++)
#add_compile_options(-stdlib=libc++)
# If it's libc++ and you're on <= 10.8, you need to compile with clang++ -stdlib=libc++. If it's libstdc++ and you're on 10.9 or later, you need to compile with clang++ -stdlib=libstdc++.
# add_compile_options(-stdlib=libstdc++)
# add_compile_options(-stdlib=libc++)
endif()
else() # if gcc
@@ -243,34 +249,14 @@ set(CMAKE_LIBRARY_PATH ${EXTERNAL_DIR}/lib $ENV{EXTERNAL_DIR}/lib ${CMAKE_LIBRAR
if (J) # dummy to avoid error msg
endif()
# set(OPTIONAL_PACKAGES "")
# list(APPEND OPTIONAL_PACKAGES "Triangle")
# list(APPEND OPTIONAL_PACKAGES "UMFPACK")
# list(APPEND OPTIONAL_PACKAGES "CHOLMOD")
# list(APPEND OPTIONAL_PACKAGES "readproc")
# list(APPEND OPTIONAL_PACKAGES "cppunit")
# list(APPEND OPTIONAL_PACKAGES "Python")
# list(APPEND OPTIONAL_PACKAGES "Sphinx")
set (UMFPACK_FOUND FALSE)
find_or_build_package(Triangle triangle LOCAL)
find_or_build_package(BLAS lapack)
find_or_build_package(LAPACK lapack)
find_or_build_package(CHOLMOD suitesparse)
#find_or_build_package(UMFPACK umfpack)
find_package(UMFPACK)
message(STATUS "************** ${UMFPACK_FOUND}")
# Check if conda package is created
if(DEFINED ENV{CONDA_BUILD})
message(STATUS "Conda package is being created.")
set(CONDA_BUILD TRUE)
else()
set(CONDA_BUILD FALSE)
endif()
if (";${BLAS_LIBRARIES};" MATCHES "openblas")
message(STATUS "openblas is used: ${BLAS_LIBRARIES}")
message(STATUS "openblas is used: ${PC_BLAS_INCLUDE_DIRS}")
@@ -449,6 +435,9 @@ endif()
message(STATUS "**********************************************************************")
message(STATUS "************************* Dependencies found *************************")
message(STATUS "**********************************************************************")
if (CONDA_BUILD)
message(STATUS "CONDA_BUILD: ${CONDA_BUILD}")
endif()
message(STATUS "THREADS :${Threads_FOUND} ${CMAKE_THREAD_LIBS_INIT}")
message(STATUS "USE_BOOST_THREAD :${USE_BOOST_THREAD} ${Boost_THREAD_LIBRARIES}")
message(STATUS "CHOLMOD_LIBRARIES :${CHOLMOD_LIBRARIES}")
@@ -459,6 +448,7 @@ message(STATUS "PYTHONLIBS_FOUND :${PYTHONLIBS_FOUND} PYTHON_LIBRARY: ${PYT
if(WIN32)
message(STATUS "PYTHON_INCLUDE_DIR :${PYTHON_INCLUDE_DIR}" )
endif(WIN32)
message(STATUS "Boost_INCLUDE_DIR :${Boost_INCLUDE_DIR}")
message(STATUS "Boost_PYTHON_FOUND :${Boost_PYTHON_FOUND} Boost_PYTHON_LIBRARY: ${Boost_PYTHON_LIBRARY}" )
message(STATUS "numpy_FOUND :${numpy_FOUND} PY_NUMPY: ${PY_NUMPY}" )
message(STATUS "CASTER_FOUND :${CASTER_FOUND} Caster: ${CASTER_EXECUTABLE}")
View
@@ -145,15 +145,24 @@ if (CHOLMOD_LIBRARY)
endif()
endif()
find_program(GFORTRAN_EXECUTABLE gfortran)
if (GFORTRAN_EXECUTABLE)
execute_process(COMMAND ${GFORTRAN_EXECUTABLE} -print-file-name=libgfortran.so
OUTPUT_VARIABLE GFORTRAN_LIBRARY
OUTPUT_STRIP_TRAILING_WHITESPACE)
if (EXISTS "${GFORTRAN_LIBRARY}")
set(CHOLMOD_LIBRARIES ${CHOLMOD_LIBRARIES} ${GFORTRAN_LIBRARY})
endif()
endif(GFORTRAN_EXECUTABLE)
if (CONDA_BUILD)
find_library(GFORTRAN_LIBRARY gfortran
DOC "The gfortran library"
)
set(CHOLMOD_LIBRARIES ${CHOLMOD_LIBRARIES} ${GFORTRAN_LIBRARY})
else()
find_program(GFORTRAN_EXECUTABLE gfortran)
if (GFORTRAN_EXECUTABLE)
execute_process(COMMAND ${GFORTRAN_EXECUTABLE} -print-file-name=libgfortran.so
OUTPUT_VARIABLE GFORTRAN_LIBRARY
OUTPUT_STRIP_TRAILING_WHITESPACE)
if (EXISTS "${GFORTRAN_LIBRARY}")
set(CHOLMOD_LIBRARIES ${CHOLMOD_LIBRARIES} ${GFORTRAN_LIBRARY})
endif()
endif(GFORTRAN_EXECUTABLE)
endif()
IF(WIN32)
View
@@ -77,6 +77,10 @@ if(WIN32)
include_directories(${PYTHON_INCLUDE_DIR})
endif(WIN32)
if(CONDA_BUILD)
include_directories(${Boost_INCLUDE_DIR})
endif(CONDA_BUILD)
ADD_CUSTOM_TARGET(pgtest)
ADD_CUSTOM_COMMAND(
@@ -483,7 +483,7 @@ def simulate(mesh, slowness, scheme, verbose=False, **kwargs):
if not ret.allNonZero('err'):
ret.set('t', t)
ret.set('err', pg.physics.Refraction.estimateError(
ret, absoluteError=noiseAbs)
ret, absoluteError=noiseAbs))
if verbose:
print("Data error estimates (min:max) ",
@@ -146,25 +146,25 @@ void DataContainerERT::checkDataValidityLocal(){
}
}
SIndex DataContainerERT::electrodeToCurrentPattern(SIndex a, SIndex b) const {
Index DataContainerERT::electrodeToCurrentPattern(Index a, Index b) const {
//std::cout << a << " " << b << " " << electrodes_.size() << std::endl;
/*! +1 ensures that -1 (pol) electrodes are considered. */
return (a + 1) * (sensorCount() + 1) + b + 1; //** better take min(a,b) and max(a,b)?
}
CurrentPattern DataContainerERT::currentPatternToElectrode(SIndex pattern){
CurrentPattern DataContainerERT::currentPatternToElectrode(Index pattern){
/*! +1 ensures that -1 (pol) electrodes are considered. */
SIndex nElecs = sensorCount() + 1;
SIndex a = (SIndex)ceil(pattern / nElecs);
SIndex b = pattern - a * nElecs;
Index nElecs = sensorCount() + 1;
Index a = (SIndex)ceil(pattern / nElecs);
Index b = pattern - a * nElecs;
a-=1; b-=1;
//std::cout << pattern << " " << a << " " << b << std::endl;
return CurrentPattern(a, b);
}
std::set < SIndex > DataContainerERT::currentPattern(bool reciprocity){
std::set < SIndex > pattern;
std::set < Index > DataContainerERT::currentPattern(bool reciprocity){
std::set < Index > pattern;
const RVector & valid = dataMap_["valid"];
const RVector & a = dataMap_["a"];
@@ -67,11 +67,11 @@ class DLLEXPORT DataContainerERT : public GIMLI::DataContainer{
virtual void checkDataValidityLocal();
CurrentPattern currentPatternToElectrode(SIndex pattern);
CurrentPattern currentPatternToElectrode(Index pattern);
SIndex electrodeToCurrentPattern(SIndex a, SIndex b) const;
Index electrodeToCurrentPattern(Index a, Index b) const;
std::set < SIndex > currentPattern(bool reciprocity=false);
std::set < Index > currentPattern(bool reciprocity=false);
/*! Merge duplicate data by averaging. Sort the DataContainerERT as well.*/
void averageDuplicateData(bool verbose=false);
@@ -312,20 +312,23 @@ void dcfemBoundaryAssembleStiffnessMatrix(CSparseMatrix & S, const Mesh & mesh,
void assembleCompleteElectrodeModel_(RSparseMatrix & S,
const std::vector < ElectrodeShape * > & elecs,
uint oldMatSize, bool lastIsReferenz){
uint oldMatSize, bool lastIsReferenz,
const RVector & contactImpedances){
RSparseMapMatrix mapS(S);
ElementMatrix < double > Se;
uint nElectrodes = elecs.size();
mapS.setRows(oldMatSize + nElectrodes);
mapS.setCols(oldMatSize + nElectrodes);
std::vector < double > vContactResistance(nElectrodes, 1.0); // Ohm
std::vector < double > vContactImpedance( nElectrodes, 1.0); // Ohm * m^2
// RVector vContactImpedance( nElectrodes, 1.0); // Ohm * m^2
bool hasImp = checkIfMapFileExistAndLoadToVector("contactImpedance.map", vContactImpedance);
// bool hasImp = checkIfMapFileExistAndLoadToVector("contactImpedance.map", vContactImpedance);
bool hasImp = true;
RVector vContactResistance(nElectrodes, 1.0); // Ohm
bool hasRes = checkIfMapFileExistAndLoadToVector("contactResistance.map", vContactResistance);
for (uint elecID = 0; elecID < nElectrodes; elecID ++){
//** some scale value, can used for contact impedance
@@ -337,7 +340,7 @@ void assembleCompleteElectrodeModel_(RSparseMatrix & S,
elecs[elecID]->setMID(mat_ID);
double contactResistance = vContactResistance[elecID];
double contactImpedance = vContactImpedance[elecID];
double contactImpedance = contactImpedances[elecID];
std::vector < MeshEntity * > electrodeEnts(elecs[elecID]->entities());
if (hasImp || hasRes){
@@ -410,13 +413,15 @@ void assembleCompleteElectrodeModel_(RSparseMatrix & S,
void assembleCompleteElectrodeModel(RSparseMatrix & S,
const std::vector < ElectrodeShape * > & elecs,
uint oldMatSize, bool lastIsReferenz){
assembleCompleteElectrodeModel_(S, elecs, oldMatSize, lastIsReferenz);
uint oldMatSize, bool lastIsReferenz,
const RVector & contactImpedances){
assembleCompleteElectrodeModel_(S, elecs, oldMatSize, lastIsReferenz, contactImpedances);
}
void assembleCompleteElectrodeModel(CSparseMatrix & S,
const std::vector < ElectrodeShape * > & elecs,
uint oldMatSize, bool lastIsReferenz){
uint oldMatSize, bool lastIsReferenz,
const RVector & contactImpedances){
THROW_TO_IMPL
}
@@ -781,6 +786,10 @@ void DCMultiElectrodeModelling::updateDataDependency_(){
if (mesh_) searchElectrodes_();
}
void DCMultiElectrodeModelling::setContactImpedances(const RVector & zi){
vContactImpedance_ = zi;
}
void DCMultiElectrodeModelling::searchElectrodes_(){
if (!mesh_){
@@ -1442,15 +1451,15 @@ void DCMultiElectrodeModelling::createCurrentPattern(std::vector < ElectrodeShap
//** this is only useful for the forward calculation since the reciprocity potentials are needed for sensitivity calculation.
if (dataContainer_){
//** reciprocity disabled
std::set < SIndex > inject(this->dataContainer().currentPattern(false));
std::set < Index > inject(this->dataContainer().currentPattern(false));
if (verbose_) std::cout << "Found " << inject.size()
<< " dipole-current pattern" << std::endl;
eA.resize(inject.size(), NULL);
eB.resize(inject.size(), NULL);
CurrentPattern cp;
uint i = 0;
for (std::set < SIndex >::iterator it = inject.begin(); it != inject.end(); it ++, i++){
for (std::set < Index >::iterator it = inject.begin(); it != inject.end(); it ++, i++){
currentPatternIdxMap_[(*it)] = i;
cp = this->dataContainer().currentPatternToElectrode((*it));
if (cp.first < nElecs && cp.second < nElecs){
@@ -1778,9 +1787,23 @@ MEMINFO
for (Index i = 0; i < passiveCEM_.size(); i ++) elecs.push_back(passiveCEM_[i]);
assembleCompleteElectrodeModel(S_, elecs, oldMatSize, lastIsReferenz_);
if (vContactImpedance_.size() == 0){
vContactImpedance_.resize(elecs.size(), 1.0); // Ohm
// RVector vContactResistance(nElectrodes, 1.0); // Ohm
// RVector vContactImpedance( nElectrodes, 1.0); // Ohm * m^2
bool hasImp = checkIfMapFileExistAndLoadToVector("contactImpedance.map",
vContactImpedance_);
if (hasImp){
if (verbose_) std::cout << "Loaded: contactImpedance.map." << std::endl;
}
}
assembleCompleteElectrodeModel(S_, elecs, oldMatSize, lastIsReferenz_,
vContactImpedance_);
potentialsCEM_.resize(nCurrentPattern, lastValidElectrode);
} // end CEM
this->assembleStiffnessMatrixDCFEMByPass(S_);
@@ -1875,12 +1898,14 @@ void DCMultiElectrodeModelling::calculateK(const std::vector < ElectrodeShape *
calculateK_(eA, eB, solutionK, kIdx);
}
void DCMultiElectrodeModelling::calculateK(const std::vector < ElectrodeShape * > & eA,
const std::vector < ElectrodeShape * > & eB,
CMatrix & solutionK, int kIdx){
calculateK_(eA, eB, solutionK, kIdx);
}
void DCSRMultiElectrodeModelling::updateMeshDependency_(){
DCMultiElectrodeModelling::updateMeshDependency_();
if (primMeshOwner_ && primMesh_) delete primMesh_;
View
@@ -46,9 +46,9 @@ DLLEXPORT double mixedBoundaryCondition(const Boundary & boundary,
double k=0.0);
DLLEXPORT void assembleCompleteElectrodeModel(RSparseMatrix & S,
const std::vector < ElectrodeShape * > & elecs, uint oldMatSize, bool lastIsReferenz);
// DLLEXPORT std::vector < ElectrodeShape * > findCEMElectrodes(const Mesh & mesh);
const std::vector < ElectrodeShape * > & elecs,
uint oldMatSize, bool lastIsReferenz,
const RVector & contactImpedances);
/*!ERT utility function for the handling of complex resistivity
* values vs. amplitude/phase data.
@@ -223,7 +223,7 @@ class DLLEXPORT DCMultiElectrodeModelling : public GIMLI::ModellingBase {
/*! Return dipole current pattern map.
* corresponds to < CurrentPattern,Idx of Potentialmatrix > */
inline const std::map < long, uint > & currentPatternIdxMap() const {
inline const std::map < Index, Index > & currentPatternIdxMap() const {
return currentPatternIdxMap_;
}
@@ -246,6 +246,8 @@ class DLLEXPORT DCMultiElectrodeModelling : public GIMLI::ModellingBase {
return electrodes_;
}
void setContactImpedances(const RVector & zi);
virtual RVector calcGeometricFactor(const DataContainerERT & data,
Index nModel=0);
@@ -325,9 +327,10 @@ class DLLEXPORT DCMultiElectrodeModelling : public GIMLI::ModellingBase {
bool buildCompleteElectrodeModel_;
bool dipoleCurrentPattern_;
std::map < long, uint > currentPatternIdxMap_;
std::map < Index, Index > currentPatternIdxMap_;
RMatrix potentialsCEM_;
RVector vContactImpedance_;
DataMap * primDataMap_;
};

0 comments on commit 11b53a8

Please sign in to comment.