4 * Created on: 7 fevrier. 2012
12 #include "MEDFileMesh.hxx"
13 #include "MEDFileField1TS.hxx"
14 #include "MEDLoader.hxx"
15 #include "MEDCouplingUMesh.hxx"
16 #include "MEDCouplingFieldDouble.hxx"
18 #include "CdmathException.hxx"
24 using namespace MEDCoupling;
28 //----------------------------------------------------------------------
29 Field::Field( EntityType typeField )
30 //----------------------------------------------------------------------
34 _numberOfComponents=0;
37 //----------------------------------------------------------------------
39 //----------------------------------------------------------------------
41 //std::cerr << "dtor Field, _field = " <<_field << std::endl;
42 //if (_field) _field->decrRef();
46 Field::Field(const std::string fieldName, EntityType type, const Mesh& mesh, int numberOfComponents, double time)
51 _numberOfComponents=numberOfComponents;
55 buildFieldMemoryStructure();
57 void Field::buildFieldMemoryStructure()
59 MEDCouplingUMesh* mu=_mesh.getMEDCouplingMesh()->buildUnstructured();
60 DataArrayDouble *array=DataArrayDouble::New();
61 if (_typeField==CELLS)
63 _field=MEDCouplingFieldDouble::New(ON_CELLS);
64 array->alloc(_mesh.getNumberOfCells(),_numberOfComponents);
66 }else if(_typeField==NODES)
68 _field=MEDCouplingFieldDouble::New(ON_NODES);
69 array->alloc(mu->getNumberOfNodes(),_numberOfComponents);
71 }else if(_typeField==FACES)
73 _field=MEDCouplingFieldDouble::New(ON_CELLS);
74 array->alloc(_mesh.getNumberOfFaces(),_numberOfComponents);
75 DataArrayIdType *desc=DataArrayIdType::New();
76 DataArrayIdType *descI=DataArrayIdType::New();
77 DataArrayIdType *revDesc=DataArrayIdType::New();
78 DataArrayIdType *revDescI=DataArrayIdType::New();
79 MEDCouplingUMesh *m3=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
87 throw CdmathException("Type of Field::Field() is not compatible");
89 _field->setName(_fieldName.c_str()) ;
90 _field->setArray(array);
91 _field->setTime(_time,0,0);
95 Field::Field( const std::string filename, EntityType type,
96 const std::string & fieldName,
97 int iteration, int order, int meshLevel)
101 _fieldName=fieldName;
103 readFieldMed(filename, type, fieldName, iteration, order);
106 Field::Field(const std::string meshFileName, EntityType type, const std::vector<double> Vconstant,
107 const std::string & fieldName, int meshLevel, double time, std::string meshName )
110 _mesh=Mesh(meshFileName + ".med", meshName, meshLevel);
112 _numberOfComponents=Vconstant.size();
114 _fieldName=fieldName;
116 buildFieldMemoryStructure();
118 int nbelem=_field->getNumberOfTuples();
119 int nbcomp=_field->getNumberOfComponents() ;
121 for (int ielem=0 ; ielem<nbelem; ielem++)
122 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
123 _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp];
125 Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std::string & fieldName, double time)
130 _numberOfComponents=Vconstant.size();
132 _fieldName=fieldName;
134 buildFieldMemoryStructure();
136 int nbelem=_field->getNumberOfTuples();
137 int nbcomp=_field->getNumberOfComponents() ;
139 for (int ielem=0 ; ielem<nbelem; ielem++)
140 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
141 _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp];
143 Field::Field(const Mesh& M, EntityType type, const vector<double> Vconstant, const std::string & fieldName, double time)
148 _numberOfComponents=Vconstant.size();
150 _fieldName=fieldName;
152 buildFieldMemoryStructure();
154 int nbelem=_field->getNumberOfTuples();
155 int nbcomp=_field->getNumberOfComponents() ;
157 for (int ielem=0 ; ielem<nbelem; ielem++)
158 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
159 _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp];
161 Field::Field( int nDim, const vector<double> Vconstant, EntityType type,
162 double xmin, double xmax, int nx, string leftSide, string rightSide,
163 double ymin, double ymax, int ny, string backSide, string frontSide,
164 double zmin, double zmax, int nz, string bottomSide, string topSide,
165 const std::string & fieldName, double time,double epsilon)
169 _numberOfComponents=Vconstant.size();
171 _fieldName=fieldName;
175 _mesh=Mesh(xmin,xmax,nx);
178 _mesh=Mesh(xmin,xmax,nx,ymin,ymax,ny);
180 _mesh=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
182 cout<<"Field::Field: Space dimension nDim should be between 1 and 3"<<endl;
183 throw CdmathException("Space dimension nDim should be between 1 and 3");
186 _mesh.setGroupAtPlan(xmax,0,epsilon,rightSide);
187 _mesh.setGroupAtPlan(xmin,0,epsilon,leftSide);
189 _mesh.setGroupAtPlan(ymax,1,epsilon,frontSide);
190 _mesh.setGroupAtPlan(ymin,1,epsilon,backSide);
193 _mesh.setGroupAtPlan(zmax,2,epsilon,topSide);
194 _mesh.setGroupAtPlan(zmin,2,epsilon,bottomSide);
198 buildFieldMemoryStructure();
200 int nbelem=_field->getNumberOfTuples();
201 int nbcomp=_field->getNumberOfComponents() ;
203 for (int ielem=0 ; ielem<nbelem; ielem++)
204 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
205 _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp];
207 Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos,
208 EntityType type, int direction, const std::string & fieldName, double time)
210 if (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
211 throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
216 _numberOfComponents=VV_Left.getNumberOfRows();
218 _fieldName=fieldName;
221 buildFieldMemoryStructure();
223 int nbelem=_field->getNumberOfTuples();
224 int nbcomp=_field->getNumberOfComponents() ;
225 double component_value;
227 for (int j = 0; j < nbelem; j++) {
229 component_value=M.getCell(j).x();
230 else if(direction==1)
231 component_value=M.getCell(j).y();
232 else if(direction==2)
233 component_value=M.getCell(j).z();
235 throw CdmathException( "Field::Field: direction should be an integer between 0 and 2");
237 for (int i=0; i< nbcomp; i++)
238 if (component_value< disc_pos )
239 _field->getArray()->getPointer()[j+i*nbcomp] = VV_Left[i];
241 _field->getArray()->getPointer()[j+i*nbcomp] = VV_Right[i];
244 Field::Field( int nDim, const vector<double> VV_Left, vector<double> VV_Right,
245 double xstep, EntityType type,
246 double xmin, double xmax, int nx, string leftSide, string rightSide,
247 double ymin, double ymax, int ny, string backSide, string frontSide,
248 double zmin, double zmax, int nz, string bottomSide, string topSide,
249 int direction, const std::string & fieldName, double time, double epsilon)
251 if (VV_Right.size()!=VV_Left.size())
252 throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
255 M=Mesh(xmin,xmax,nx);
257 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
259 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
261 throw CdmathException("Field::Field : Space dimension nDim should be between 1 and 3");
263 M.setGroupAtPlan(xmax,0,epsilon,rightSide);
264 M.setGroupAtPlan(xmin,0,epsilon,leftSide);
266 M.setGroupAtPlan(ymax,1,epsilon,frontSide);
267 M.setGroupAtPlan(ymin,1,epsilon,backSide);
270 M.setGroupAtPlan(zmax,2,epsilon,topSide);
271 M.setGroupAtPlan(zmin,2,epsilon,bottomSide);
273 Vector V_Left(VV_Left.size()), V_Right(VV_Right.size());
274 for(int i=0;i<VV_Left.size(); i++){
275 V_Left(i)=VV_Left[i];
276 V_Right(i)=VV_Right[i];
278 Field(M, V_Left, V_Right, xstep, type, direction, fieldName, time);
281 Field::Field(const Mesh M, const Vector Vin, const Vector Vout, double radius,
282 const Vector Center, EntityType type, const std::string & fieldName, double time)
284 if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
286 cout<< "Vout.size()= "<<Vout.size() << ", Vin.size()= "<<Vin.size()<<", Center.size()="<<Center.size()<<", M.getSpaceDim= "<< M.getSpaceDimension()<<endl;
287 throw CdmathException("Field::Field : Vector size error");
293 _numberOfComponents=Vout.size();
295 _fieldName=fieldName;
298 buildFieldMemoryStructure();
300 int nbelem=_field->getNumberOfTuples();
301 int nbcomp=_field->getNumberOfComponents() ;
303 int spaceDim=M.getSpaceDimension();
304 Vector currentPoint(spaceDim);
306 for(int i=0;i<nbelem;i++)
308 currentPoint(0)=M.getCell(i).x();
311 currentPoint(1)=M.getCell(i).y();
313 currentPoint(2)=M.getCell(i).z();
315 if((currentPoint-Center).norm()<radius)
316 for(int j=0;j<nbcomp;j++)
317 _field->getArray()->getPointer()[j+i*nbcomp]=Vin[j];
319 for(int j=0;j<nbcomp;j++)
320 _field->getArray()->getPointer()[j+i*nbcomp]=Vout[j];
324 MEDCoupling::DataArrayDouble * Field::getArray(){
325 return _field->getArray();
329 Field::readFieldMed( const std::string & fileNameRadical,
331 const std::string & fieldName,
336 * Reads the file fileNameRadical.med and creates a Field from it.
338 std::string completeFileName = fileNameRadical + ".med";
339 std::vector<std::string> fieldNames = MEDCoupling::GetAllFieldNames(completeFileName);
341 std::string attributedFieldName;
344 // Get the name of the right field that we will attribute to the Field.
345 if (fieldName == "") {
346 if (fieldNames.size() > 0)
348 cout<<"Warning : No field name imposed, taking the first field name found : "<< fieldNames[0]<<" in file "<<completeFileName<<endl;;
349 attributedFieldName = fieldNames[0];
353 std::ostringstream message;
354 message << "Warning : No field found in file " << completeFileName;
355 throw CdmathException(message.str().c_str());
360 int mycount = std::count(fieldNames.begin(), fieldNames.end(), fieldName);
364 attributedFieldName = fieldName;
367 cout<<"Warning : " << mycount << " fields are associated to the name " << fieldName <<" in file "<<completeFileName<<endl;
370 std::ostringstream message;
371 message << "No field named " << fieldName << " in file " << completeFileName;
372 throw CdmathException(message.str().c_str());
376 // Get the name of the right mesh that we will attribute to the Field.
377 std::vector<std::string> meshNames = MEDCoupling::GetMeshNamesOnField(completeFileName, attributedFieldName);
378 if (meshNames.size() == 0) {
379 std::ostringstream message;
380 message << "No mesh associated to " << fieldName<< " in file " << completeFileName;
381 throw CdmathException(message.str().c_str());
384 if( meshNames.size() > 1 )
386 cout<<"Warning : " << meshNames.size() << " meshes are associated to field named " << fieldName <<" in file "<<completeFileName<<endl;
387 cout<<"Mesh names are : ";
388 for (int iMesh=0; iMesh < meshNames.size(); iMesh++)
389 cout<< meshNames[iMesh]<<", ";
390 cout<<endl<<"Taking the mesh with name "<< meshNames[0]<<endl;
393 std::string attributedMeshName = meshNames[0];
399 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
400 MEDCoupling::ReadFieldCell( completeFileName,
401 attributedMeshName, 0,
402 attributedFieldName, iteration, order) );
405 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
406 MEDCoupling::ReadFieldNode( completeFileName,
407 attributedMeshName, 0,
408 attributedFieldName, iteration, order) );
411 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
412 MEDCoupling::ReadFieldCell( completeFileName,
413 attributedMeshName, -1,
414 attributedFieldName, iteration, order) );
418 //Read and store the number of components
419 _numberOfComponents = _field->getNumberOfComponents() ;
420 _time = _field->getTime(iteration, order);
422 cout<<"Found field " << fieldName << " in file " << completeFileName << " on mesh named "<< _field->getMesh()->getName()<< endl;
423 _mesh=Mesh( completeFileName, _field->getMesh()->getName());
428 Field::getNormEuclidean() const
430 Vector result(_numberOfComponents);
431 DoubleTab norm(_numberOfComponents,_field->magnitude()->getArray()->getConstPointer());
432 result.setValues(norm);
438 Field::max(int component) const
440 if( component >= getNumberOfComponents() || component < 0)
441 throw CdmathException("double Field::max() : component number should be between 0 and the field number of components");
443 double result=-1e100;
444 for(int i=0; i<getNumberOfElements() ; i++)
445 if( result < (*this)(i,component))
446 result = (*this)(i,component);
452 Field::min(int component) const
454 if( component >= getNumberOfComponents() || component < 0)
455 throw CdmathException("double Field::min() : component number should be between 0 and the field number of components");
458 for(int i=0; i<getNumberOfElements() ; i++)
459 if( result > (*this)(i,component))
460 result = (*this)(i,component);
466 Field::getInfoOnComponent(int icomp) const
468 return _field->getArray()->getInfoOnComponent(icomp);
472 Field::setInfoOnComponent(int icomp, string nameCompo)
474 _field.retn()->getArray()->setInfoOnComponent(icomp,nameCompo);
477 //----------------------------------------------------------------------
478 Field::Field( const Field & f )
479 //----------------------------------------------------------------------
482 MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
484 _typeField=f.getTypeOfField();
487 //----------------------------------------------------------------------
488 MCAuto<MEDCouplingFieldDouble>
489 Field::getField ( void ) const
490 //----------------------------------------------------------------------
495 //----------------------------------------------------------------------
497 Field::setFieldByMEDCouplingFieldDouble ( const MEDCouplingFieldDouble* field )
498 //----------------------------------------------------------------------
500 MCAuto<MEDCouplingFieldDouble> ff=field->deepCopy();
504 //----------------------------------------------------------------------
506 Field::setFieldByDataArrayDouble ( const DataArrayDouble* array )
507 //----------------------------------------------------------------------
509 _field->setArray(const_cast<DataArrayDouble*>(array));
512 //----------------------------------------------------------------------
514 Field::integral ( ) const
515 //----------------------------------------------------------------------
517 int nbComp=_field->getNumberOfComponents();
518 double res[nbComp];//Pointer containing the inegral of each component
519 _field->integral(false,res);
520 Vector result(nbComp);//Vector containing the inegral of each component
522 for(int i=0; i<nbComp ; i++)
528 //----------------------------------------------------------------------
530 Field::integral ( int compId ) const
531 //----------------------------------------------------------------------
533 return _field->integral(compId, false);
536 //----------------------------------------------------------------------
538 Field::normL1 ( ) const
539 //----------------------------------------------------------------------
541 int nbComp=_field->getNumberOfComponents();
542 double res[nbComp];//Pointer containing the L1 norm of each component
544 Vector result(nbComp);//Vector containing the L1 norm of each component
546 for(int i=0; i<nbComp ; i++)
552 //----------------------------------------------------------------------
554 Field::normL2 ( ) const
555 //----------------------------------------------------------------------
557 int nbComp=_field->getNumberOfComponents();
558 double res[nbComp];//Pointer containing the L2 norm of each component
560 Vector result(nbComp);//Vector containing the L2 norm of each component
562 for(int i=0; i<nbComp ; i++)
568 //----------------------------------------------------------------------
570 Field::normMax ( ) const
571 //----------------------------------------------------------------------
573 int nbComp=_field->getNumberOfComponents();
574 double res[nbComp];//Pointer containing the L2 norm of each component
575 _field->normMax(res);
576 Vector result(nbComp);//Vector containing the L2 norm of each component
578 for(int i=0; i<nbComp ; i++)
585 Field::componentMax(Vector & Indices) const
587 int nbComp=_field->getNumberOfComponents();
588 int nbElems=getNumberOfElements();
590 Vector result(nbComp);//Vector containing the Linfinity norm of each component
592 for(int i=0; i<nbElems ; i++)
593 for(int j=0; j<nbComp ; j++)
594 if(fabs((*this)(i,j))>result(j))
596 result(j)=fabs((*this)(i,j));
602 //----------------------------------------------------------------------
604 Field::operator() ( int ielem )
605 //----------------------------------------------------------------------
607 if(ielem>_field->getNumberOfTuples() || ielem<0)
608 throw CdmathException("double& Field::operator(ielem) : ielem>number of values !");
609 return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
612 //----------------------------------------------------------------------
614 Field::operator[] ( int ielem )
615 //----------------------------------------------------------------------
617 if(ielem>_field->getNumberOfTuples() || ielem<0)
618 throw CdmathException("double& Field::operator[ielem] : ielem>number of values !");
619 return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
622 //----------------------------------------------------------------------
624 Field::operator() ( int ielem ) const
625 //----------------------------------------------------------------------
627 if(ielem>_field->getNumberOfTuples() || ielem<0)
628 throw CdmathException("double Field::operator(ielem) : ielem>number of values !");
629 return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
632 //----------------------------------------------------------------------
634 Field::operator[] ( int ielem ) const
635 //----------------------------------------------------------------------
637 if(ielem>_field->getNumberOfTuples() || ielem<0)
638 throw CdmathException("double Field::operator[ielem] : ielem>number of values !");
639 return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
642 //----------------------------------------------------------------------
644 Field::operator() ( int ielem, int jcomp )
645 //----------------------------------------------------------------------
647 if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
648 throw CdmathException("double& Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
649 return _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()];
652 //----------------------------------------------------------------------
654 Field::operator() ( int ielem, int jcomp ) const
655 //----------------------------------------------------------------------
657 if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
658 throw CdmathException("double Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
659 return _field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()];
662 //----------------------------------------------------------------------
664 Field::setTime ( double time, int iter )
665 //----------------------------------------------------------------------
667 _field->setTime(time,iter,0.0);
669 //----------------------------------------------------------------------
671 Field::getTime ( void ) const
672 //----------------------------------------------------------------------
675 return _field->getTime(a,b);
678 //----------------------------------------------------------------------
680 Field::getNumberOfElements ( void ) const
681 //----------------------------------------------------------------------
683 return _field->getNumberOfTuples() ;
687 Field::getSpaceDimension( void ) const
689 return _mesh.getSpaceDimension() ;
692 //----------------------------------------------------------------------
694 Field::getNumberOfComponents ( void ) const
695 //----------------------------------------------------------------------
697 return _field->getNumberOfComponents() ;
700 //----------------------------------------------------------------------
702 Field::getValues ( void ) const
703 //----------------------------------------------------------------------
705 return _field->getArray()->getConstPointer() ;
708 //----------------------------------------------------------------------
710 Field::getValues ( Vector myVector ) const
711 //----------------------------------------------------------------------
713 if( myVector.size() != _field->getNumberOfTuples() * _field->getNumberOfComponents() )
714 throw CdmathException("Error : Field::getValues requires a vector having the same number of component as fiedl values");
716 const double * fieldValues = _field->getArray()->getConstPointer();
717 double * vectorValues = myVector.getValues().getValues();
719 memcpy(vectorValues, fieldValues,myVector.size()*sizeof(double)) ;
723 Field::setValues ( Vector myVector )
724 //----------------------------------------------------------------------
726 if( myVector.size() != _field->getNumberOfTuples() * _field->getNumberOfComponents() )
727 throw CdmathException("Error : Field::setValues requires a vector having the same number of component as fiedl values");
729 double * vectorValues = myVector.getValues().getValues();
731 setValues ( vectorValues, myVector.size() );
735 Field::setValues ( double * values, int nbValues )
737 double * fieldValues = _field->getArray()->getPointer() ;
738 memcpy(fieldValues,values,nbValues*sizeof(double)) ;
741 //----------------------------------------------------------------------
743 Field::getName ( void ) const
744 //----------------------------------------------------------------------
746 return _field->getName() ;
749 //----------------------------------------------------------------------
751 Field::getMesh ( void ) const
752 //----------------------------------------------------------------------
757 //----------------------------------------------------------------------
759 Field::getTypeOfField ( void ) const
760 //----------------------------------------------------------------------
766 Field::getElementComponent(int i, int comp) const
774 return _mesh.getCell(i).x();
776 return _mesh.getCell(i).y();
778 return _mesh.getCell(i).z();
780 cout<<"Wrong component number "<< comp <<" , dimension is "<< _mesh.getSpaceDimension() << ", field values are on CELLS" <<endl;
781 throw CdmathException("Field::getElementComponent : Wrong component number");
787 return _mesh.getNode(i).x();
789 return _mesh.getNode(i).y();
791 return _mesh.getNode(i).z();
793 cout<<"Wrong component number "<< comp <<" , dimension is "<< _mesh.getSpaceDimension() << ", field values are on NODES" <<endl;
794 throw CdmathException("Field::getElementComponent : Wrong component number");
800 return _mesh.getFace(i).x();
802 return _mesh.getFace(i).y();
804 return _mesh.getFace(i).z();
806 cout<<"Wrong component number "<< comp <<" , dimension is "<< _mesh.getSpaceDimension() << ", field values are on FACES"<<endl;
807 throw CdmathException("Field::getElementComponent : Wrong component number");
810 throw CdmathException("Accepted field supports are CELLS, NODES and FACES");
814 //----------------------------------------------------------------------
816 Field::setName ( const string fieldName )
817 //----------------------------------------------------------------------
819 _field->setName(fieldName.c_str()) ;
823 //----------------------------------------------------------------------
825 Field::operator+ ( const Field& f ) const
826 //----------------------------------------------------------------------
828 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
829 //throw CdmathException("Field::operator+ : Field addition requires identical meshes");
830 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
831 if(f.getTypeOfField() != getTypeOfField())
832 throw CdmathException("Field::operator+ : Field addition requires identical field types (CELLS, NODES or FACES");
833 if(f.getNumberOfComponents() != getNumberOfComponents())
834 throw CdmathException("Field::operator+ : Field addition requires identical number of components");
836 Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
837 int nbComp=f.getNumberOfComponents();
838 int nbElem=f.getNumberOfElements();
839 for (int ielem=0 ; ielem<nbElem; ielem++)
840 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
841 fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*nbComp]+f(ielem, jcomp);
845 //----------------------------------------------------------------------
847 Field::operator- ( const Field& f ) const
848 //----------------------------------------------------------------------
850 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
851 //throw CdmathException("Field::operator- : Field subtraction requires identical meshes");
852 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
853 if(f.getTypeOfField() != getTypeOfField())
854 throw CdmathException("Field::operator- : Field subtraction requires identical field types (CELLS, NODES or FACES");
855 if(f.getNumberOfComponents() != getNumberOfComponents())
856 throw CdmathException("Field::operator- : Field subtraction requires identical number of components");
858 Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
859 int nbComp=f.getNumberOfComponents();
860 int nbElem=f.getNumberOfElements();
861 for (int ielem=0 ; ielem<nbElem; ielem++)
862 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
863 fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*nbComp]-f(ielem, jcomp);
867 //----------------------------------------------------------------------
869 Field::operator= ( const Field& f )
870 //----------------------------------------------------------------------
873 _typeField=f.getTypeOfField() ;
874 _numberOfComponents=f.getNumberOfComponents();
876 _fieldName=f.getName();
877 MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
882 Field Field::deepCopy( ) const
884 Field F(getName(), getTypeOfField(), getMesh(), getNumberOfComponents(), getTime()) ;
885 MCAuto<MEDCouplingFieldDouble> f1=getField()->deepCopy();
886 F.setFieldByMEDCouplingFieldDouble(f1);
892 //----------------------------------------------------------------------
894 Field::operator+= ( const Field& f )
895 //----------------------------------------------------------------------
897 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
898 //throw CdmathException("Field::operator+= : Field addition requires identical meshes");
899 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
900 if(f.getTypeOfField() != getTypeOfField())
901 throw CdmathException("Field::operator+= : Field addition requires identical field types (CELLS, NODES or FACES");
902 if(f.getNumberOfComponents() != getNumberOfComponents())
903 throw CdmathException("Field::operator+= : Field addition requires identical number of components");
905 _field->setMesh(f.getField()->getMesh());
906 (*_field)+=(*f.getField());
910 //----------------------------------------------------------------------
912 Field::operator-= ( const Field& f )
913 //----------------------------------------------------------------------
915 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
916 //throw CdmathException("Field::operator-= : Field subtraction requires identical meshes");
917 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
918 if(f.getTypeOfField() != getTypeOfField())
919 throw CdmathException("Field::operator-= : Field subtraction requires identical field types (CELLS, NODES or FACES");
920 if(f.getNumberOfComponents() != getNumberOfComponents())
921 throw CdmathException("Field::operator-= : Field subtraction requires identical number of components");
923 _field->setMesh(f.getField()->getMesh());
924 (*_field)-=(*f.getField());
928 //----------------------------------------------------------------------
930 Field::operator*= ( double s )
931 //----------------------------------------------------------------------
933 int nbComp=getNumberOfComponents();
934 int nbElem=getNumberOfElements();
935 for (int i=0 ; i<nbComp ; i++)
936 for (int j=0 ; j<nbElem; j++)
937 _field->getArray()->getPointer()[i+j*nbComp]*=s;
941 //----------------------------------------------------------------------
943 Field::operator/= ( double s )
944 //----------------------------------------------------------------------
946 int nbComp=getNumberOfComponents();
947 int nbElem=getNumberOfElements();
948 for (int i=0 ; i<nbComp ; i++)
949 for (int j=0 ; j<nbElem; j++)
950 _field->getArray()->getPointer()[i+j*nbComp]/=s;
954 //----------------------------------------------------------------------
956 Field::operator-= ( double s )
957 //----------------------------------------------------------------------
959 int nbComp=getNumberOfComponents();
960 int nbElem=getNumberOfElements();
961 for (int i=0 ; i<nbComp ; i++)
962 for (int j=0 ; j<nbElem; j++)
963 _field->getArray()->getPointer()[i+j*nbComp]-=s;
967 //----------------------------------------------------------------------
969 Field::operator+= ( double s )
970 //----------------------------------------------------------------------
972 int nbComp=getNumberOfComponents();
973 int nbElem=getNumberOfElements();
974 for (int i=0 ; i<nbComp ; i++)
975 for (int j=0 ; j<nbElem; j++)
976 _field->getArray()->getPointer()[i+j*nbComp]+=s;
980 //----------------------------------------------------------------------
982 Field::writeVTK (std::string fileName, bool fromScratch) const
983 //----------------------------------------------------------------------
985 string fname=fileName+".pvd";
987 double time=_field->getTime(iter,order);
991 ofstream file(fname.c_str()) ;
992 file << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\"><Collection>\n" ;
993 ostringstream numfile;
995 string filetmp=fileName+"_";
996 filetmp=filetmp+numfile.str();
997 string ret=_field->writeVTK(filetmp.c_str()) ;
998 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
999 file << "</Collection></VTKFile>\n" ;
1004 ifstream file1(fname.c_str()) ;
1006 getline(file1, contenus, '\0');
1007 string to_remove="</Collection></VTKFile>";
1008 size_t m = contenus.find(to_remove);
1009 size_t n = contenus.find_first_of("\n", m + to_remove.length());
1010 contenus.erase(m, n - m + 1);
1012 ofstream file(fname.c_str()) ;
1014 ostringstream numfile;
1016 string filetmp=fileName+"_";
1017 filetmp=filetmp+numfile.str();
1018 string ret=_field->writeVTK(filetmp.c_str()) ;
1019 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
1020 file << "</Collection></VTKFile>\n" ;
1025 //----------------------------------------------------------------------
1027 Field::writeCSV ( const std::string fileName ) const
1028 //----------------------------------------------------------------------
1031 double time=_field->getTime(iter,order);
1033 ostringstream numfile;
1035 string filetmp=fileName+"_";
1036 filetmp=filetmp+numfile.str();
1037 filetmp=filetmp+".csv";
1038 ofstream file(filetmp.c_str()) ;
1039 int dim=_mesh.getSpaceDimension();
1041 if (getTypeOfField()==CELLS)
1042 nbElements=_mesh.getNumberOfCells();
1044 nbElements=_mesh.getNumberOfNodes();
1048 int nbCompo=getNumberOfComponents();
1050 file << "x " << _field->getName() << endl;
1054 for (int i=0;i<nbCompo;i++)
1055 file << " " << _field->getName() << " Compo " << i+1 << " "<< getInfoOnComponent(i);
1058 for (int i=0;i<nbElements;i++)
1060 if (getTypeOfField()==CELLS)
1061 file << _mesh.getCell(i).x() ;
1063 file << _mesh.getNode(i).x() ;
1064 for (int j=0;j<nbCompo;j++)
1065 file << " " << getValues()[j+i*nbCompo] ;
1070 int nbCompo=getNumberOfComponents();
1072 file << "x y " << _field->getName() << endl;
1076 for (int i=0;i<nbCompo;i++)
1077 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
1080 for (int i=0;i<nbElements;i++)
1082 if (getTypeOfField()==CELLS)
1083 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() ;
1085 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() ;
1086 for (int j=0;j<nbCompo;j++)
1087 file << " " << getValues()[j+i*nbCompo] ;
1092 int nbCompo=getNumberOfComponents();
1094 file << "x y z " << _field->getName() << endl;
1098 for (int i=0;i<nbCompo;i++)
1099 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
1102 for (int i=0;i<nbElements;i++)
1104 if (getTypeOfField()==CELLS)
1105 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() << " " << _mesh.getCell(i).z();
1107 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() << " " << _mesh.getNode(i).z();
1108 for (int j=0;j<nbCompo;j++)
1109 file << " " << getValues()[j+i*nbCompo] ;
1116 //----------------------------------------------------------------------
1118 Field::writeMED ( const std::string fileName, bool fromScratch)
1119 //----------------------------------------------------------------------
1121 string fname=fileName+".med";
1123 if(_mesh.isStructured() || _mesh.meshNotDeleted())
1125 MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
1127 MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
1128 else//The mesh has ben deleted, use a MEDFileField1TS instead of _field to save the values
1130 if ( not fromScratch)
1132 MEDCoupling::MCAuto<MEDCoupling::MEDFileField1TS> ff=MEDFileField1TS::New();//To save the field when the mesh has been deleted
1133 ff->setFieldNoProfileSBT( _field );
1134 ff->write(fname.c_str(), fromScratch);
1137 throw CdmathException("Field::writeMED Error !!! The mesh has been deleted, cannot write field from scratch");
1142 operator* (double value , const Field& field )
1144 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1145 int nbComp=field.getNumberOfComponents();
1146 int nbElem=field.getNumberOfElements();
1147 for (int ielem=0 ; ielem<nbElem; ielem++)
1148 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1149 fres(ielem, jcomp)=value*field(ielem, jcomp);
1154 operator* (const Field& field, double value )
1156 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1157 int nbComp=field.getNumberOfComponents();
1158 int nbElem=field.getNumberOfElements();
1159 for (int ielem=0 ; ielem<nbElem; ielem++)
1160 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1161 fres(ielem, jcomp)=value*field(ielem, jcomp);
1165 Field operator/ (const Field& field, double value)
1167 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1168 int nbComp=field.getNumberOfComponents();
1169 int nbElem=field.getNumberOfElements();
1170 for (int ielem=0 ; ielem<nbElem; ielem++)
1171 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1172 fres(ielem, jcomp)=field(ielem, jcomp)/value;
1177 Field::getValuesOnAllComponents(int elem) const
1179 Vector v(getNumberOfComponents());
1180 for(int i=0;i<getNumberOfComponents();i++)
1181 v[i]=(*this)(elem,i);
1186 Field::getValuesOnComponent(int compo) const
1188 Vector v(getNumberOfElements());
1189 for(int i=0;i<getNumberOfElements();i++)
1190 v[i]=(*this)(i,compo);
1194 std::vector< double >
1195 Field::getFieldValues(int compo) const
1197 std::vector< double > v(getNumberOfElements());
1198 for(int i=0;i<getNumberOfElements();i++)
1199 v[i]=(*this)(i,compo);
1203 std::ostream& operator<<(std::ostream& out, const Field& field )
1205 cout << "Field " << field.getName() << " : " << endl ;
1206 out<< field.getField().retn()->getArray()->repr();
1210 void Field::deleteMEDCouplingMesh()
1212 return _mesh.deleteMEDCouplingMesh();