recvcount, recvtype, source, recvtag, comm,status);
}
int worldSize() const {int size; MPI_Comm_size(MPI_COMM_WORLD, &size); return size;}
+ int reduce(void* sendbuf, void* recvbuf, int count, MPI_Datatype datatype, MPI_Op op,
+ int root, MPI_Comm comm) const
+ {return MPI_Reduce(sendbuf, recvbuf, count, datatype, op, root, comm);}
+ int allReduce(void* sendbuf, void* recvbuf, int count, MPI_Datatype datatype, MPI_Op op,
+ MPI_Comm comm) const
+ {return MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm);}
};
}
#include "ComponentTopology.hxx"
#include "ParaFIELD.hxx"
#include "DEC.hxx"
+#include "./ICoCo/ICoCoField.hxx"
+#include "./ICoCo/ICoCoMEDField.hxx"
+#include "./ICoCo/ICoCoTrioField.hxx"
namespace ParaMEDMEM
{
DEC::DEC(ProcessorGroup& source_group, ProcessorGroup& target_group):_local_field(0),
- _source_group(&source_group), _target_group(&target_group)
+ _source_group(&source_group), _target_group(&target_group), _owns_field(false)
{
_union_group = source_group.fuse(target_group);
}
DEC::~DEC()
{
delete _union_group;
+ if (_owns_field)
+ delete _local_field;
}
-void DEC::attachLocalField(const ParaFIELD* field)
-{
- _local_field=field;
- //if (field!=0)
- //{
- //BlockTopology* topo=dynamic_cast<BlockTopology*>(field->getTopology());
- _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
- //}
-}
+ void DEC::attachLocalField(const ParaFIELD* field)
+ {
+ _local_field=field;
+ _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
+ }
+
+ void DEC::attachLocalField(const ICoCo::Field* field)
+ {
+ const ICoCo::MEDField* medfield =dynamic_cast<const ICoCo::MEDField*> (field);
+ if(medfield !=0)
+ {
+ attachLocalField(medfield->getField());
+ return;
+ }
+ const ICoCo::TrioField* triofield=dynamic_cast<const ICoCo::TrioField*> (field);
+ if (triofield !=0)
+ {
+ ProcessorGroup* localgroup;
+ if (_source_group->containsMyRank()) localgroup=_source_group;
+ else localgroup=_target_group;
+
+ ICoCo::Field* parafield=new ICoCo::MEDField(*triofield, *localgroup);
+ _owns_field = true;
+ attachLocalField(parafield);
+ }
+ throw MEDMEM::MEDEXCEPTION("incompatible field type");
+ }
+ void DEC::setForcedRenormalizationFlag(bool flag)
+ {
+ _forced_renormalization_flag=flag;
+ }
+ bool DEC::getForcedRenormalizationFlag()
+ {
+ return _forced_renormalization_flag;
+ }
}
#ifndef DEC_HXX_
#define DEC_HXX_
-
+namespace ICoCo
+{
+ class Field;
+}
namespace ParaMEDMEM
{
class ProcessorGroup;
DEC():_local_field(0){}
DEC(ProcessorGroup& local_group, ProcessorGroup& distant_group);
void attachLocalField(const ParaFIELD* field);
+ void attachLocalField(const ICoCo::Field* field);
virtual void prepareSourceDE()=0;
virtual void prepareTargetDE()=0;
virtual void recvData()=0;
virtual void synchronize()=0;
virtual ~DEC();
virtual void computeProcGroup(){};
+ void setForcedRenormalizationFlag(bool flag);
+ bool getForcedRenormalizationFlag();
protected:
const ParaFIELD* _local_field;
//! Processor group representing the union of target and source processors
ProcessorGroup* _target_group;
const CommInterface* _comm_interface;
+ bool _forced_renormalization_flag;
+ bool _owns_field;
};
}
*/
void InterpolationMatrix::addContribution(MEDMEM::MESH& distant_support, int iproc_distant, int* distant_elems)
{
+
+ if (distant_support.getMeshDimension() != _source_support.getMeshDimension() ||
+ distant_support.getMeshDimension() != _source_support.getMeshDimension() )
+ throw MEDMEM::MEDEXCEPTION("local and distant meshes do not have the same space and mesh dimensions");
+
//creating the interpolator structure
- MEDMEM::Interpolation* interpolator;
+ MEDMEM::Interpolation* interpolator;
if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==3)
interpolator=new MEDMEM::Interpolation3DSurf();
else if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==2)
interpolator=new MEDMEM::Interpolation2D();
+ else
+ throw MEDMEM::MEDEXCEPTION("no interpolator exists for these mesh and space dimensions ");
//computation of the intersection volumes between source and target elements
int source_size= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
throw MEDEXCEPTION("uncoherent number of rows in interpolation matrix");
//computing the vectors containing the source and target element volumes
- MEDMEM::SUPPORT target_support(&distant_support,"all cells", MED_EN::MED_CELL);
- MEDMEM::FIELD<double>* target_triangle_surf = distant_support.getArea(&target_support);
+ MEDMEM::SUPPORT target_support(&distant_support,"all cells", MED_EN::MED_CELL);
+ MEDMEM::FIELD<double>* target_triangle_surf = getSupportVolumes(target_support);
MEDMEM::SUPPORT source_support (const_cast<MEDMEM::MESH*>(&_source_support),"all cells", MED_EN::MED_CELL);
- MEDMEM::FIELD<double>* source_triangle_surf = _source_support.getArea(&source_support);
+ MEDMEM::FIELD<double>* source_triangle_surf = getSupportVolumes(source_support);
//storing the source volumes
_source_volume.resize(source_size);
//VS^(-1) is the inverse of the diagonal matrix storing
//volumes of the source cells
for (int i=0; i<nbrows; i++)
- for (int icomp; icomp<nbcomp; icomp++)
+ for (int icomp=0; icomp<nbcomp; icomp++)
target_value[i*nbcomp+icomp]/=_source_volume[i];
//storing S= VS^(-1).(WT.T)
}
+/*!
+\brief returns the volumes of the cells underlying the field \a field
+For 2D geometries, the returned field contains the areas.
+For 3D geometries, the returned field contains the volumes.
+
+\param field field on which cells the volumes are required
+\return field containing the volumes
+*/
+ MEDMEM::FIELD<double>* InterpolationMatrix::getSupportVolumes(const MEDMEM::SUPPORT& support)
+ {
+ const MEDMEM::MESH* mesh=support.getMesh();
+ int dim = mesh->getMeshDimension();
+ switch (dim)
+ {
+ case 2:
+ return mesh->getArea(&support);
+ case 3:
+ return mesh->getVolume(&support);
+ default:
+ throw MEDMEM::MEDEXCEPTION("interpolation is not available for this dimension");
+ }
+ }
}
const ProcessorGroup& distant_group,
const string& method);
- //InterpolationMatrix(const MEDMEM::MESH& source_support, const string& method);
- //InterpolationMatrix(const MEDMEM::MESH& target_support, const MEDMEM::MESH& source_support);
virtual ~InterpolationMatrix();
void addContribution(MEDMEM::MESH& distant_support, int iproc_distant, int* distant_elems);
void multiply(MEDMEM::FIELD<double>&) const;
int getNbRows() const {return _row_offsets.size();}
private:
- // vector<pair <int,int> > _source_indices;
+
+ MEDMEM::FIELD<double>* InterpolationMatrix::getSupportVolumes(const MEDMEM::SUPPORT& field);
+
+ private:
vector<int> _row_offsets;
- //vector<int> _col_numbers;
vector<pair<int,int> > _col_offsets;
- // vector<double> _coeffs;
const MEDMEM::MESH& _source_support;
MxN_Mapping _mapping;
string _method;
}
}
_interpolation_matrix->prepare();
- // _volume_vector=para_mesh->getVolume();
}
if (_source_group->containsMyRank())
_interpolation_matrix->transposeMultiply(*_local_field->getField());
else if (_target_group->containsMyRank())
- _interpolation_matrix->multiply(*_local_field->getField());
+ {
+ _interpolation_matrix->multiply(*_local_field->getField());
+ if (_forced_renormalization_flag)
+ for (int icomp=0; icomp<_local_field->getField()->getNumberOfComponents(); icomp++)
+ {
+ double total_norm = _local_field->getVolumeIntegral(icomp+1);
+ double source_norm=total_norm;
+ _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
+
+ if (_target_group->containsMyRank() && abs(total_norm)>1e-100)
+ _local_field->getField()->applyLin(source_norm/total_norm,0.0,icomp+1);
+ }
+ }
+
+
}
void IntersectionDEC::sendData()
{
if (_source_group->containsMyRank())
- _interpolation_matrix->multiply(*_local_field->getField());
+ {
+ _interpolation_matrix->multiply(*_local_field->getField());
+ if (_forced_renormalization_flag)
+ for (int icomp=0; icomp<_local_field->getField()->getNumberOfComponents(); icomp++)
+ {
+ double total_norm = _local_field->getVolumeIntegral(icomp+1);
+ double source_norm = total_norm;
+ _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
+
+ if (_target_group->containsMyRank() && abs(total_norm)>1e-100)
+ _local_field->getField()->applyLin(source_norm/total_norm,0.0,icomp+1);
+ }
+ }
else if (_target_group->containsMyRank())
_interpolation_matrix->transposeMultiply(*_local_field->getField());
}
top_srcdir=@top_srcdir@
top_builddir=../..
srcdir=@srcdir@
+
VPATH=.:$(srcdir):$(srcdir)/tests
MACHINE=PCLINUX
@COMMENCE@
+SUBDIRS= ICoCo
EXPORT_PYSCRIPTS = \
#include "UnstructuredParaSUPPORT.hxx"
#include "ExplicitCoincidentDEC.hxx"
#include "StructuredCoincidentDEC.hxx"
-
+#include "CommInterface.hxx"
+#include "ProcessorGroup.hxx"
+#include "MPIProcessorGroup.hxx"
#include "ParaFIELD.hxx"
#include "ParaMESH.hxx"
}
}
-// int nb_components=0;
-// if (component_topology.getProcGroup()!=0)
+
int nb_components = component_topology.nbLocalComponents();
if (nb_components!=0)
{
}
else return;
- _field->setName("toto");
- _field->setDescription("titi");
+ _field->setName("Default ParaFIELD name");
+ _field->setDescription("Default ParaFIELD description");
_field->setNumberOfValues(para_support->getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS));
string* compnames=new string[nb_components];
string* compdesc=new string[nb_components];
delete data_channel;
}
+double ParaFIELD::getVolumeIntegral(int icomp) const
+{
+ CommInterface comm_interface = _topology->getProcGroup()->getCommInterface();
+ const MEDMEM::SUPPORT* support = _field->getSupport();
+ FIELD<double>* volume= getSupportVolumes(*support);
+ double integral=0;
+ for (int i=0; i<support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS); i++)
+ integral+=_field->getValueIJ(i+1,icomp)*volume->getValueIJ(i+1,1);
+
+ double total=0.0;
+ comm_interface.allReduce(&integral, &total, 1, MPI_DOUBLE, MPI_SUM, * dynamic_cast<const MPIProcessorGroup*>(_topology->getProcGroup())->getComm());
+ return total;
+}
+
+/*!
+\brief returns the volumes of the cells underlying the field \a field
+
+For 2D geometries, the returned field contains the areas.
+For 3D geometries, the returned field contains the volumes.
+
+\param field field on which cells the volumes are required
+\return field containing the volumes
+*/
+ MEDMEM::FIELD<double>* ParaFIELD::getSupportVolumes(const MEDMEM::SUPPORT& support) const
+ {
+ const MEDMEM::MESH* mesh=support.getMesh();
+ int dim = mesh->getMeshDimension();
+ switch (dim)
+ {
+ case 2:
+ return mesh->getArea(&support);
+ case 3:
+ return mesh->getVolume(&support);
+ default:
+ throw MEDMEM::MEDEXCEPTION("interpolation is not available for this dimension");
+ }
+ }
}
const ParaSUPPORT* getSupport() const {return _support;}
Topology* getTopology() const {return _topology;}
int nbComponents() const {return _component_topology.nbComponents();}
+ double getVolumeIntegral(int icomp) const;
+ double getL2norm()const{return -1;}
private:
+ MEDMEM::FIELD<double>* getSupportVolumes(const MEDMEM::SUPPORT& support) const;
const ComponentTopology& _component_topology;
Topology* _topology;
MEDMEM::FIELD<double>* _field;
{
}
- ParaSUPPORT::ParaSUPPORT(const MEDMEM::SUPPORT& support, const ProcessorGroup& proc_group):_support(&support) {
+ ParaSUPPORT::ParaSUPPORT(const MEDMEM::SUPPORT& support, const ProcessorGroup& proc_group):_support(&support), _owns_support(false) {
_mesh = new ParaMESH(*(support.getMesh()), proc_group, "mesh from support");
}
ParaSUPPORT::~ParaSUPPORT()
{
+ if (_owns_support)
+ delete _support;
}
}
public:
ParaSUPPORT();
ParaSUPPORT(const ParaMESH* mesh, const MEDMEM::SUPPORT* support):
- _support(support), _mesh(mesh){};
- ParaSUPPORT(const MEDMEM::SUPPORT&, const ProcessorGroup& proc_group);
- virtual ~ParaSUPPORT();
+ _support(support), _mesh(mesh), _owns_support(false){};
+ ParaSUPPORT(const ParaMESH* mesh);
+ ParaSUPPORT(const MEDMEM::SUPPORT&, const ProcessorGroup& proc_group);
+ virtual ~ParaSUPPORT();
virtual const Topology* getTopology() const{};
virtual const MEDMEM::SUPPORT* getSupport() const {return _support;}
virtual const ParaMESH* getMesh() const {return _mesh;}
protected :
+
const MEDMEM::SUPPORT* _support;
const ParaMESH* _mesh;
+ bool _owns_support;
};
}
-#include "INTERPOLATION_2D.hxx"
+#include "Interpolation2D.hxx"
#include "MEDMEM_Mesh.hxx"
int main()
MEDMEM::MESH source(MED_DRIVER,"/home/vb144235/resources/square128000.med","Mesh_1");
MEDMEM::MESH target(MED_DRIVER,"/home/vb144235/resources/square30000.med","Mesh_1");
- ParaMEDMEM::INTERPOLATION_2D interp;
+ MEDMEM::Interpolation2D interp;
// interp.setOptions(1e-6,1,2);
interp.interpol_maillages(source,target);
}
//loading the geometry for the source group
if (source_group->containsMyRank())
{
- string master = "/home/vb144235/resources/square128000_split";
+ // string master = "/home/vb144235/resources/square128000_split";
+ string master = "/home/vb144235/resources/square1_split";
+
cout <<"loading source"<<endl;
source_mesh=new ParaMESH(MED_DRIVER,master,*source_group);
cout <<"end of load"<<endl;
strstream <<master<<rank+1<<".med";
//strstream<<"/home/vb144235/resources/square128000.med";
ostringstream meshname ;
- meshname<< "Mesh_1_"<< rank+1;
+ meshname<< "Mesh_2_"<< rank+1;
// meshname<<"Mesh_1";
MEDMEM::MESH* mesh_bis = new MESH(MED_DRIVER,strstream.str(),meshname.str());
source_mesh_bis = new ParaMESH(*mesh_bis,*source_group,"mesh_from_mem");
//loading the geometry for the target group
if (target_group->containsMyRank())
{
- string master = "/home/vb144235/resources/square30000_split";
+ // string master = "/home/vb144235/resources/square30000_split";
+ string master= "/home/vb144235/resources/square2_split";
target_mesh=new ParaMESH(MED_DRIVER,master,*target_group);
topo_target=target_mesh->getBlockTopology();
ostringstream strstream;
strstream << master<<(rank-nproc_source+1)<<".med";
//strstream<<"/home/vb144235/resources/square30000.med";
ostringstream meshname ;
- meshname<< "Mesh_1_"<<rank-nproc_source+1;
+ meshname<< "Mesh_3_"<<rank-nproc_source+1;
//meshname<<"Mesh_1";
MEDMEM::MESH* mesh_bis = new MESH(MED_DRIVER,strstream.str(),meshname.str());
source_field.getField()->setValue(value);
cout <<"creating intersectionDEC"<<endl;
IntersectionDEC dec(*source_group,*target_group);
- dec.attachSourceField(&source_field);
+ dec.attachLocalField(&source_field);
dec.synchronize();
cout<<"DEC usage"<<endl;
+ dec.setForcedRenormalizationFlag(true);
+
dec.sendData();
+ source_mesh->write(MED_DRIVER,"/home/vb144235/tmp/sourcesquareb");
+ source_field.write(MED_DRIVER,"/home/vb144235/tmp/sourcesquareb","boundary");
+
+ dec.recvData();
cout <<"writing"<<endl;
source_mesh->write(MED_DRIVER,"/home/vb144235/tmp/sourcesquare");
source_field.write(MED_DRIVER,"/home/vb144235/tmp/sourcesquare","boundary");
value[ielem]=0.0;
target_field.getField()->setValue(value);
IntersectionDEC dec(*source_group,*target_group);
- dec.attachTargetField(&target_field);
+ dec.attachLocalField(&target_field);
dec.synchronize();
+ dec.setForcedRenormalizationFlag(true);
+
dec.recvData();
+ target_mesh->write(MED_DRIVER, "/home/vb144235/tmp/targetsquareb");
+ target_field.write(MED_DRIVER, "/home/vb144235/tmp/targetsquareb", "boundary");
+ dec.sendData();
target_mesh->write(MED_DRIVER, "/home/vb144235/tmp/targetsquare");
target_field.write(MED_DRIVER, "/home/vb144235/tmp/targetsquare", "boundary");
delete[] value;
source_field.getField()->setValue(value);
cout <<"creating intersectionDEC"<<endl;
IntersectionDEC dec(*source_group,*target_group);
- dec.attachSourceField(&source_field);
+ dec.attachLocalField(&source_field);
dec.synchronize();
cout<<"DEC usage"<<endl;
dec.sendData();
value[ielem]=0.0;
target_field.getField()->setValue(value);
IntersectionDEC dec(*source_group,*target_group);
- dec.attachTargetField(&target_field);
+ dec.attachLocalField(&target_field);
dec.synchronize();
dec.recvData();
target_mesh->write(MED_DRIVER, "/home/vb144235/tmp/targetsquare");
source_field.getField()->setValue(value);
cout <<"creating intersectionDEC"<<endl;
IntersectionDEC dec(*source_group,*target_group);
- dec.attachSourceField(&source_field);
+ dec.attachLocalField(&source_field);
dec.synchronize();
cout<<"DEC usage"<<endl;
dec.sendData();
value[ielem]=0.0;
target_field.getField()->setValue(value);
IntersectionDEC dec(*source_group,*target_group);
- dec.attachTargetField(&target_field);
+ dec.attachLocalField(&target_field);
dec.synchronize();
dec.recvData();
target_mesh_bis->write(MED_DRIVER, "/home/vb144235/tmp/targetcylinder");
source_field.getField()->setValue(value);
cout <<"creating intersectionDEC"<<endl;
IntersectionDEC dec(*source_group,*target_group);
- dec.attachSourceField(&source_field);
+ dec.attachLocalField(&source_field);
dec.synchronize();
cout<<"DEC usage"<<endl;
dec.sendData();
+ dec.recvData();
cout <<"writing"<<endl;
source_mesh->write(MED_DRIVER,"/home/vb144235/tmp/sourcesquare");
source_field.write(MED_DRIVER,"/home/vb144235/tmp/sourcesquare","boundary");
value[ielem]=0.0;
target_field.getField()->setValue(value);
IntersectionDEC dec(*source_group,*target_group);
- dec.attachTargetField(&target_field);
+ dec.attachLocalField(&target_field);
dec.synchronize();
dec.recvData();
+ dec.sendData();
target_mesh->write(MED_DRIVER, "/home/vb144235/tmp/targetsquare");
target_field.write(MED_DRIVER, "/home/vb144235/tmp/targetsquare", "boundary");
delete[] value;