#include "Interpolation3D.txx"
#include "MEDNormalizedUnstructuredMesh.hxx"
#include "InterpolationOptions.hxx"
+#include "IntersectionResult.hxx"
/*! \class InterpolationMatrix
_row_offsets[i]=0;
_coeffs.resize(nbelems);
+ _coef_values.reserve( target_group.size() );
}
InterpolationMatrix::~InterpolationMatrix()
int iproc_distant, int* distant_elems)
{
if (distant_support.getMeshDimension() != _source_support.getMeshDimension() ||
- distant_support.getMeshDimension() != _source_support.getMeshDimension() )
+ distant_support.getSpaceDimension() != _source_support.getSpaceDimension() )
throw MEDMEM::MEDEXCEPTION("local and distant meshes do not have the same space and mesh dimensions");
+ bool isP1Method = ( getMethod() == "P1" || getMethod() == "P1d" );
+
//creating the interpolator structure
- vector<map<int,double> > surfaces;
+ //vector<map<int,double> > surfaces;
+ INTERP_KERNEL::IntersectionResult surfaces( _source_support.getSpaceDimension() );
+ surfaces.setNeedBaryCentres( isP1Method );
+
//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);
MEDNormalizedUnstructuredMesh<2,2> source_wrapper(&_source_support);
INTERP_KERNEL::Interpolation2D interpolator (*this);
+ if ( distant_support.getNumberOfPolygons(MED_CELL) > 0 ||
+ _source_support.getNumberOfPolygons(MED_CELL) > 0 )
+ // it must be DualMESH
+ interpolator.setIntersectionType( INTERP_KERNEL::Geometric2D );
interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces);
}
else if (distant_support.getMeshDimension()==3 && distant_support.getSpaceDimension()==3)
"all cells", MED_EN::MED_CELL);
MEDMEM::FIELD<double>* source_triangle_surf = getSupportVolumes(source_support);
+ // allocate vector of coefficients
+ int nbNodesPerCell = source_support.getTypes()[0] % 100;
+ int nbValuesPerCell = isP1Method ? nbNodesPerCell : 1;
+ int nbIntersections = surfaces.getNumberOfIntersections();
+ _coef_values.push_back(vector<double>());
+ vector<double>& coef_values = _coef_values.back();
+ coef_values.reserve( nbIntersections * nbValuesPerCell );
+
+ // prepare to getting coordinates of nodes of a source cell
+ // needed for comuting barycentric coordinates of intersection gravity center
+ const int* src_conn = _source_support.getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,
+ MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ const double * src_coords = _source_support.getCoordinates(MED_EN::MED_FULL_INTERLACE);
+ vector<const double*> cell_coords( nbNodesPerCell );
+ vector<double> bary_coords( nbNodesPerCell );
+
//storing the source volumes
_source_volume.resize(source_size);
for (int i=0; i<source_size; i++)
//col_id is the number of the column
//iter->second is the value of the coefficient
- _coeffs[ielem].push_back(make_pair(col_id,iter->second));
+ //_coeffs[ielem].push_back(make_pair(col_id,iter->second));
+ if ( nbValuesPerCell == 1 )
+ {
+ coef_values.push_back(iter->second);
+ _coeffs[ielem].push_back(make_pair(col_id, &coef_values.back() ));
+ }
+ else
+ {
+ // barycenter of intersection of cells ielem and iter->first
+ double* bary_centre = surfaces.getBaryCentre( ielem, iter->first );
+ // coordinates of a source cell
+ const int* cell_conn = src_conn + nbNodesPerCell * ielem;
+ for ( int i = 0; i < nbNodesPerCell; ++i, ++cell_conn )
+ cell_coords[i] = src_coords + _source_support.getSpaceDimension() * ( *cell_conn-1 );
+ // barycentric coodrs
+ INTERP_KERNEL::barycentric_coodrs( cell_coords, bary_centre, &bary_coords[0]);
+ // store (surface * bary_coords)
+ coef_values.push_back(iter->second * bary_coords[0]);
+ _coeffs[ielem].push_back(make_pair(col_id, &coef_values.back() ));
+ for ( int j = 1; j < bary_coords.size(); ++j )
+ coef_values.push_back(iter->second * bary_coords[j]);
+ }
}
}
delete source_triangle_surf;
delete target_triangle_surf;
}
-
+
/*! The call to this method updates the arrays on the target side
so that they know which amount of data from which processor they
should expect.
// W is the intersection matrix
// S is the source vector
- for (int irow=0; irow<nbrows; irow++)
- for (int icomp=0; icomp< nbcomp; icomp++)
- {
- double coeff_row = field.getValueIJ(irow+1,icomp+1);
- for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+ if ( getMethod() == "P0" ) { // --------------- P0 - one value per cell
+
+ for (int irow=0; irow<nbrows; irow++)
+ for (int icomp=0; icomp< nbcomp; icomp++)
{
- int colid= _coeffs[irow][icol-_row_offsets[irow]].first;
- double value = _coeffs[irow][icol-_row_offsets[irow]].second;
- target_value[(colid-1)*nbcomp+icomp]+=value*coeff_row;
+ double coeff_row = field.getValueIJ(irow+1,icomp+1);
+ for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+ {
+ int colid = _coeffs[irow][icol-_row_offsets[irow]].first;
+ double* value = _coeffs[irow][icol-_row_offsets[irow]].second;
+ target_value[(colid-1)*nbcomp+icomp]+=value[0]*coeff_row;
+ }
}
- }
-
+ }
+ else if ( getMethod() == "P1" ) { // ---------- P1 - one value per node
+
+ int nbValuesPerCell = _source_support.getTypes(MED_EN::MED_CELL)[0] % 100;
+ const int* cell_conn = _source_support.getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,
+ MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+
+ for (int irow=0; irow<nbrows; irow++, cell_conn += nbValuesPerCell)
+ for (int icomp=0; icomp< nbcomp; icomp++)
+ for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+ {
+ int colid = _coeffs[irow][icol-_row_offsets[irow]].first;
+ double* value = _coeffs[irow][icol-_row_offsets[irow]].second;
+ for ( int ival=0; ival<nbValuesPerCell; ++ival)
+ {
+ double coeff_row = field.getValueIJ(cell_conn[ival],icomp+1);
+ target_value[(colid-1)*nbcomp+icomp]+=value[ival]*coeff_row;
+ }
+ }
+ }
+ else { // ------------------------------------- P1d - different values per node
+
+ int nbValuesPerCell = _source_support.getTypes(MED_EN::MED_CELL)[0] % 100;
+
+ for (int irow=0; irow<nbrows; irow++)
+ for (int icomp=0; icomp< nbcomp; icomp++)
+ for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+ {
+ int colid = _coeffs[irow][icol-_row_offsets[irow]].first;
+ double* value = _coeffs[irow][icol-_row_offsets[irow]].second;
+ for ( int ival=0; ival<nbValuesPerCell; ++ival)
+ {
+ double coeff_row = field.getValueIJK(irow+1,icomp+1,ival+1);
+ target_value[(colid-1)*nbcomp+icomp]+=value[ival]*coeff_row;
+ }
+ }
+ }
+
// performing VT^(-1).(W.S)
// where VT^(-1) is the inverse of the diagonal matrix containing
// the volumes of target cells
for (int icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
{
int colid = _coeffs[irow][icol-_row_offsets[irow]].first;
- double value = _coeffs[irow][icol-_row_offsets[irow]].second;
+ double* value = _coeffs[irow][icol-_row_offsets[irow]].second;
for (int icomp=0; icomp<nbcomp; icomp++)
{
double coeff_row = source_value[(colid-1)*nbcomp+icomp];
- target_value[irow*nbcomp+icomp] += value*coeff_row;
+ target_value[irow*nbcomp+icomp] += value[0]*coeff_row;
}
}
for (int icomp=0; icomp<nbcomp; icomp++)
target_value[i*nbcomp+icomp] /= _source_volume[i];
- //storing S= VS^(-1).(WT.T)
- for (int irow=0; irow<nbrows; irow++)
- for (int icomp=0; icomp<nbcomp; icomp++)
- field.setValueIJ(irow+1,icomp+1,target_value[irow*nbcomp+icomp]);
+ if ( getMethod() == "P0" ) {
+
+ //storing S= VS^(-1).(WT.T)
+ for (int irow=0; irow<nbrows; irow++)
+ for (int icomp=0; icomp<nbcomp; icomp++)
+ field.setValueIJ(irow+1,icomp+1,target_value[irow*nbcomp+icomp]);
+
+ }
}
}
