4 * Created on: 7 fevrier. 2012
11 #include "MEDFileMesh.hxx"
12 #include "MEDLoader.hxx"
13 #include "MEDCouplingUMesh.hxx"
14 #include "MEDCouplingFieldDouble.hxx"
16 #include "CdmathException.hxx"
20 using namespace MEDCoupling;
24 //----------------------------------------------------------------------
25 Field::Field( EntityType typeField )
26 //----------------------------------------------------------------------
30 _numberOfComponents=0;
33 //----------------------------------------------------------------------
35 //----------------------------------------------------------------------
37 //std::cerr << "dtor Field, _field = " <<_field << std::endl;
38 //if (_field) _field->decrRef();
42 Field::Field(const std::string fieldName, EntityType type, const Mesh& mesh, int numberOfComponents, double time)
47 _numberOfComponents=numberOfComponents;
51 buildFieldMemoryStructure();
53 void Field::buildFieldMemoryStructure()
55 MEDCouplingUMesh* mu=_mesh.getMEDCouplingMesh()->buildUnstructured();
56 DataArrayDouble *array=DataArrayDouble::New();
57 if (_typeField==CELLS)
59 _field=MEDCouplingFieldDouble::New(ON_CELLS);
60 array->alloc(_mesh.getNumberOfCells(),_numberOfComponents);
62 }else if(_typeField==NODES)
64 _field=MEDCouplingFieldDouble::New(ON_NODES);
65 array->alloc(_mesh.getNumberOfNodes(),_numberOfComponents);
67 }else if(_typeField==FACES)
69 _field=MEDCouplingFieldDouble::New(ON_CELLS);
70 array->alloc(_mesh.getNumberOfFaces(),_numberOfComponents);
71 DataArrayInt *desc=DataArrayInt::New();
72 DataArrayInt *descI=DataArrayInt::New();
73 DataArrayInt *revDesc=DataArrayInt::New();
74 DataArrayInt *revDescI=DataArrayInt::New();
75 MEDCouplingUMesh *m3=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
83 throw CdmathException("Type of Field::Field() is not compatible");
85 _field->setName(_fieldName.c_str()) ;
86 _field->setArray(array);
87 _field->setTime(_time,0,0);
92 Field::Field( const std::string filename, EntityType type,
93 const std::string & fieldName,
94 int iteration, int order, int meshLevel,
95 int numberOfComponents, double time)
98 _mesh=Mesh(filename + ".med", meshLevel);
100 _numberOfComponents=numberOfComponents;
102 _fieldName=fieldName;
104 readFieldMed(filename, type, fieldName, iteration, order);
107 Field::Field(const std::string meshFileName, EntityType type, const std::vector<double> Vconstant,
108 const std::string & fieldName, int meshLevel, double time )
111 _mesh=Mesh(meshFileName + ".med", meshLevel);
113 _numberOfComponents=Vconstant.size();
115 _fieldName=fieldName;
117 buildFieldMemoryStructure();
119 int nbelem=_field->getNumberOfTuples();
120 int nbcomp=_field->getNumberOfComponents() ;
122 for (int ielem=0 ; ielem<nbelem; ielem++)
123 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
124 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
126 Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std::string & fieldName, double time)
131 _numberOfComponents=Vconstant.size();
133 _fieldName=fieldName;
135 buildFieldMemoryStructure();
137 int nbelem=_field->getNumberOfTuples();
138 int nbcomp=_field->getNumberOfComponents() ;
140 for (int ielem=0 ; ielem<nbelem; ielem++)
141 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
142 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
144 Field::Field(const Mesh& M, EntityType type, const vector<double> Vconstant, const std::string & fieldName, double time)
149 _numberOfComponents=Vconstant.size();
151 _fieldName=fieldName;
153 buildFieldMemoryStructure();
155 int nbelem=_field->getNumberOfTuples();
156 int nbcomp=_field->getNumberOfComponents() ;
158 for (int ielem=0 ; ielem<nbelem; ielem++)
159 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
160 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
162 Field::Field( int nDim, const vector<double> Vconstant, EntityType type,
163 double xmin, double xmax, int nx, string leftSide, string rightSide,
164 double ymin, double ymax, int ny, string backSide, string frontSide,
165 double zmin, double zmax, int nz, string bottomSide, string topSide,
166 const std::string & fieldName, double time,double epsilon)
170 _numberOfComponents=Vconstant.size();
172 _fieldName=fieldName;
176 _mesh=Mesh(xmin,xmax,nx);
179 _mesh=Mesh(xmin,xmax,nx,ymin,ymax,ny);
181 _mesh=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
183 cout<<"Field::Field: Space dimension nDim should be between 1 and 3"<<endl;
184 throw CdmathException("Space dimension nDim should be between 1 and 3");
187 _mesh.setGroupAtPlan(xmax,0,epsilon,rightSide);
188 _mesh.setGroupAtPlan(xmin,0,epsilon,leftSide);
190 _mesh.setGroupAtPlan(ymax,1,epsilon,frontSide);
191 _mesh.setGroupAtPlan(ymin,1,epsilon,backSide);
194 _mesh.setGroupAtPlan(zmax,2,epsilon,topSide);
195 _mesh.setGroupAtPlan(zmin,2,epsilon,bottomSide);
199 buildFieldMemoryStructure();
201 int nbelem=_field->getNumberOfTuples();
202 int nbcomp=_field->getNumberOfComponents() ;
204 for (int ielem=0 ; ielem<nbelem; ielem++)
205 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
206 _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
208 Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos,
209 EntityType type, int direction, const std::string & fieldName, double time)
211 if (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
212 throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
217 _numberOfComponents=VV_Left.getNumberOfRows();
219 _fieldName=fieldName;
222 buildFieldMemoryStructure();
224 int nbelem=_field->getNumberOfTuples();
225 int nbcomp=_field->getNumberOfComponents() ;
226 double component_value;
228 for (int j = 0; j < nbelem; j++) {
230 component_value=M.getCell(j).x();
231 else if(direction==1)
232 component_value=M.getCell(j).y();
233 else if(direction==2)
234 component_value=M.getCell(j).z();
236 throw CdmathException( "Field::Field: direction should be an integer between 0 and 2");
238 for (int i=0; i< nbcomp; i++)
239 if (component_value< disc_pos )
240 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Left[i];
242 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Right[i];
245 Field::Field( int nDim, const vector<double> VV_Left, vector<double> VV_Right,
246 double xstep, EntityType type,
247 double xmin, double xmax, int nx, string leftSide, string rightSide,
248 double ymin, double ymax, int ny, string backSide, string frontSide,
249 double zmin, double zmax, int nz, string bottomSide, string topSide,
250 int direction, const std::string & fieldName, double time, double epsilon)
252 if (VV_Right.size()!=VV_Left.size())
253 throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
256 M=Mesh(xmin,xmax,nx);
258 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
260 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
262 throw CdmathException("Field::Field : Space dimension nDim should be between 1 and 3");
264 M.setGroupAtPlan(xmax,0,epsilon,rightSide);
265 M.setGroupAtPlan(xmin,0,epsilon,leftSide);
267 M.setGroupAtPlan(ymax,1,epsilon,frontSide);
268 M.setGroupAtPlan(ymin,1,epsilon,backSide);
271 M.setGroupAtPlan(zmax,2,epsilon,topSide);
272 M.setGroupAtPlan(zmin,2,epsilon,bottomSide);
274 Vector V_Left(VV_Left.size()), V_Right(VV_Right.size());
275 for(int i=0;i<VV_Left.size(); i++){
276 V_Left(i)=VV_Left[i];
277 V_Right(i)=VV_Right[i];
279 Field(M, V_Left, V_Right, xstep, type, direction, fieldName, time);
282 Field::Field(const Mesh M, const Vector Vin, const Vector Vout, double radius,
283 const Vector Center, EntityType type, const std::string & fieldName, double time)
285 if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
287 cout<< "Vout.size()= "<<Vout.size() << ", Vin.size()= "<<Vin.size()<<", Center.size()="<<Center.size()<<", M.getSpaceDim= "<< M.getSpaceDimension()<<endl;
288 throw CdmathException("Field::Field : Vector size error");
294 _numberOfComponents=Vout.size();
296 _fieldName=fieldName;
299 buildFieldMemoryStructure();
301 int nbelem=_field->getNumberOfTuples();
302 int nbcomp=_field->getNumberOfComponents() ;
304 int spaceDim=M.getSpaceDimension();
305 Vector currentPoint(spaceDim);
307 for(int i=0;i<nbelem;i++)
309 currentPoint(0)=M.getCell(i).x();
312 currentPoint(1)=M.getCell(i).y();
314 currentPoint(2)=M.getCell(i).z();
316 if((currentPoint-Center).norm()<radius)
317 for(int j=0;j<nbcomp;j++)
318 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vin[j];
320 for(int j=0;j<nbcomp;j++)
321 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vout[j];
325 MEDCoupling::DataArrayDouble * Field::getArray(){
326 return _field->getArray();
330 Field::readFieldMed( const std::string & fileNameRadical,
332 const std::string & fieldName,
337 * Reads the file fileNameRadical.med and creates a Field from it.
339 std::string completeFileName = fileNameRadical + ".med";
340 std::vector<std::string> fieldNames = MEDCoupling::GetAllFieldNames(completeFileName);
342 std::string attributedFieldName;
345 // Get the name of the right field that we will attribute to the Field.
346 if (fieldName == "") {
347 if (fieldNames.size() > 0)
348 attributedFieldName = fieldNames[0];
350 std::ostringstream message;
351 message << "No field in file " << completeFileName;
352 throw CdmathException(message.str().c_str());
356 for (; iField < fieldNames.size(); iField++)
357 if (fieldName == fieldNames[iField]) break;
359 if (iField < fieldNames.size())
360 attributedFieldName = fieldName;
362 std::ostringstream message;
363 message << "No field named " << fieldName << " in file " << completeFileName;
364 throw CdmathException(message.str().c_str());
368 // Get the name of the right mesh that we will attribute to the Field.
369 std::vector<std::string> meshNames
370 = MEDCoupling::GetMeshNamesOnField(completeFileName, attributedFieldName);
371 if (meshNames.size() == 0) {
372 std::ostringstream message;
373 message << "No mesh associated to " << fieldName
374 << " in file " << completeFileName;
375 throw CdmathException(message.str().c_str());
377 std::string attributedMeshName = meshNames[0];
380 MEDCoupling::TypeOfField medFieldType[3] = { ON_CELLS, ON_NODES, ON_CELLS };
383 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
384 MEDCoupling::ReadFieldCell( completeFileName,
385 attributedMeshName, 0,
386 attributedFieldName, iteration, order) );
389 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
390 MEDCoupling::ReadFieldNode( completeFileName,
391 attributedMeshName, 0,
392 attributedFieldName, iteration, order) );
395 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
396 MEDCoupling::ReadFieldCell( completeFileName,
397 attributedMeshName, -1,
398 attributedFieldName, iteration, order) );
405 Field::getNormEuclidean() const
407 DoubleTab norm(getNumberOfElements(),_field->magnitude()->getArray()->getConstPointer());
412 Field::max(int component) const
414 if( component >= getNumberOfComponents() )
415 throw CdmathException("double Field::max() : component number should be smaller than field number of components");
417 double result=-1e100;
418 for(int i=0; i<getNumberOfElements() ; i++)
419 if( result < (*this)(i,component))
420 result = (*this)(i,component);
426 Field::min(int component) const
428 if( component >= getNumberOfComponents() )
429 throw CdmathException("double Field::min() : component number should be smaller than field number of components");
432 for(int i=0; i<getNumberOfElements() ; i++)
433 if( result > (*this)(i,component))
434 result = (*this)(i,component);
440 Field::getInfoOnComponent(int icomp) const
442 return _field->getArray()->getInfoOnComponent(icomp);
446 Field::setInfoOnComponent(int icomp, string nameCompo)
448 _field.retn()->getArray()->setInfoOnComponent(icomp,nameCompo);
451 //----------------------------------------------------------------------
452 Field::Field( const Field & f )
453 //----------------------------------------------------------------------
456 MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
458 _typeField=f.getTypeOfField();
461 //----------------------------------------------------------------------
462 MCAuto<MEDCouplingFieldDouble>
463 Field::getField ( void ) const
464 //----------------------------------------------------------------------
469 //----------------------------------------------------------------------
471 Field::setFieldByMEDCouplingFieldDouble ( const MEDCouplingFieldDouble* field )
472 //----------------------------------------------------------------------
474 MCAuto<MEDCouplingFieldDouble> ff=field->deepCopy();
478 //----------------------------------------------------------------------
480 Field::setFieldByDataArrayDouble ( const DataArrayDouble* array )
481 //----------------------------------------------------------------------
483 _field->setArray(const_cast<DataArrayDouble*>(array));
486 //----------------------------------------------------------------------
488 Field::integral ( ) const
489 //----------------------------------------------------------------------
491 int nbComp=_field->getNumberOfComponents();
492 double res[nbComp];//Pointer containing the inegral of each component
493 _field->integral(false,res);
494 Vector result(nbComp);//Vector containing the inegral of each component
496 for(int i=0; i<nbComp ; i++)
502 //----------------------------------------------------------------------
504 Field::integral ( int compId ) const
505 //----------------------------------------------------------------------
507 return _field->integral(compId, false);
510 //----------------------------------------------------------------------
512 Field::normL1 ( ) const
513 //----------------------------------------------------------------------
515 int nbComp=_field->getNumberOfComponents();
516 double res[nbComp];//Pointer containing the L1 norm of each component
518 Vector result(nbComp);//Vector containing the L1 norm of each component
520 for(int i=0; i<nbComp ; i++)
526 //----------------------------------------------------------------------
528 Field::normL2 ( ) const
529 //----------------------------------------------------------------------
531 int nbComp=_field->getNumberOfComponents();
532 double res[nbComp];//Pointer containing the L2 norm of each component
534 Vector result(nbComp);//Vector containing the L2 norm of each component
536 for(int i=0; i<nbComp ; i++)
542 //----------------------------------------------------------------------
544 Field::normMax ( ) const
545 //----------------------------------------------------------------------
547 int nbComp=_field->getNumberOfComponents();
548 int nbElems=getNumberOfElements();
550 Vector result(nbComp);//Vector containing the Linfinity norm of each component
552 for(int i=0; i<nbElems ; i++)
553 for(int j=0; j<nbComp ; j++)
554 if(fabs((*this)(i,j))>result(j))
555 result(j)=fabs((*this)(i,j));
561 Field::componentMax(Vector & Indices) const
563 int nbComp=_field->getNumberOfComponents();
564 int nbElems=getNumberOfElements();
566 Vector result(nbComp);//Vector containing the Linfinity norm of each component
568 for(int i=0; i<nbElems ; i++)
569 for(int j=0; j<nbComp ; j++)
570 if(fabs((*this)(i,j))>result(j))
572 result(j)=fabs((*this)(i,j));
578 //----------------------------------------------------------------------
580 Field::operator() ( int ielem )
581 //----------------------------------------------------------------------
583 if(ielem>_field->getNumberOfTuples() || ielem<0)
584 throw CdmathException("double& Field::operator(ielem) : ielem>number of values !");
585 return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
588 //----------------------------------------------------------------------
590 Field::operator[] ( int ielem )
591 //----------------------------------------------------------------------
593 if(ielem>_field->getNumberOfTuples() || ielem<0)
594 throw CdmathException("double& Field::operator[ielem] : ielem>number of values !");
595 return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
598 //----------------------------------------------------------------------
600 Field::operator() ( int ielem ) const
601 //----------------------------------------------------------------------
603 if(ielem>_field->getNumberOfTuples() || ielem<0)
604 throw CdmathException("double Field::operator(ielem) : ielem>number of values !");
605 return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
608 //----------------------------------------------------------------------
610 Field::operator[] ( int ielem ) const
611 //----------------------------------------------------------------------
613 if(ielem>_field->getNumberOfTuples() || ielem<0)
614 throw CdmathException("double Field::operator[ielem] : ielem>number of values !");
615 return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
618 //----------------------------------------------------------------------
620 Field::operator() ( int ielem, int jcomp )
621 //----------------------------------------------------------------------
623 if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
624 throw CdmathException("double& Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
625 return _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()];
628 //----------------------------------------------------------------------
630 Field::operator() ( int ielem, int jcomp ) const
631 //----------------------------------------------------------------------
633 if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
634 throw CdmathException("double Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
635 return _field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()];
638 //----------------------------------------------------------------------
640 Field::setTime ( double time, int iter )
641 //----------------------------------------------------------------------
643 _field->setTime(time,iter,0.0);
645 //----------------------------------------------------------------------
647 Field::getTime ( void ) const
648 //----------------------------------------------------------------------
651 return _field->getTime(a,b);
654 //----------------------------------------------------------------------
656 Field::getNumberOfElements ( void ) const
657 //----------------------------------------------------------------------
659 return _field->getNumberOfTuples() ;
663 Field::getSpaceDimension( void ) const
665 return _mesh.getSpaceDimension() ;
668 //----------------------------------------------------------------------
670 Field::getNumberOfComponents ( void ) const
671 //----------------------------------------------------------------------
673 return _field->getNumberOfComponents() ;
676 //----------------------------------------------------------------------
678 Field::getValues ( void ) const
679 //----------------------------------------------------------------------
681 return _field->getArray()->getConstPointer() ;
684 //----------------------------------------------------------------------
686 Field::getName ( void ) const
687 //----------------------------------------------------------------------
689 return _field->getName() ;
692 //----------------------------------------------------------------------
694 Field::getMesh ( void ) const
695 //----------------------------------------------------------------------
700 //----------------------------------------------------------------------
702 Field::getTypeOfField ( void ) const
703 //----------------------------------------------------------------------
708 //----------------------------------------------------------------------
710 Field::setName ( const string fieldName )
711 //----------------------------------------------------------------------
713 _field->setName(fieldName.c_str()) ;
717 //----------------------------------------------------------------------
719 Field::operator+ ( const Field& f ) const
720 //----------------------------------------------------------------------
722 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
723 //throw CdmathException("Field::operator+ : Field addition requires identical meshes");
724 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
725 if(f.getTypeOfField() != getTypeOfField())
726 throw CdmathException("Field::operator+ : Field addition requires identical field types (CELLS, NODES or FACES");
727 if(f.getNumberOfComponents() != getNumberOfComponents())
728 throw CdmathException("Field::operator+ : Field addition requires identical number of components");
730 Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
731 int nbComp=f.getNumberOfComponents();
732 int nbElem=f.getNumberOfElements();
733 for (int ielem=0 ; ielem<nbElem; ielem++)
734 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
735 fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]+f(ielem, jcomp);
739 //----------------------------------------------------------------------
741 Field::operator- ( const Field& f ) const
742 //----------------------------------------------------------------------
744 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
745 //throw CdmathException("Field::operator- : Field subtraction requires identical meshes");
746 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
747 if(f.getTypeOfField() != getTypeOfField())
748 throw CdmathException("Field::operator- : Field subtraction requires identical field types (CELLS, NODES or FACES");
749 if(f.getNumberOfComponents() != getNumberOfComponents())
750 throw CdmathException("Field::operator- : Field subtraction requires identical number of components");
752 Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
753 int nbComp=f.getNumberOfComponents();
754 int nbElem=f.getNumberOfElements();
755 for (int ielem=0 ; ielem<nbElem; ielem++)
756 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
757 fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]-f(ielem, jcomp);
761 //----------------------------------------------------------------------
763 Field::operator= ( const Field& f )
764 //----------------------------------------------------------------------
767 _typeField=f.getTypeOfField() ;
768 _numberOfComponents=f.getNumberOfComponents();
770 _fieldName=f.getName();
771 MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
776 Field Field::deepCopy( ) const
778 Field F(getName(), getTypeOfField(), getMesh(), getNumberOfComponents(), getTime()) ;
779 MCAuto<MEDCouplingFieldDouble> f1=getField()->deepCopy();
780 F.setFieldByMEDCouplingFieldDouble(f1);
786 //----------------------------------------------------------------------
788 Field::operator+= ( const Field& f )
789 //----------------------------------------------------------------------
791 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
792 //throw CdmathException("Field::operator+= : Field addition requires identical meshes");
793 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
794 if(f.getTypeOfField() != getTypeOfField())
795 throw CdmathException("Field::operator+= : Field addition requires identical field types (CELLS, NODES or FACES");
796 if(f.getNumberOfComponents() != getNumberOfComponents())
797 throw CdmathException("Field::operator+= : Field addition requires identical number of components");
799 _field->setMesh(f.getField()->getMesh());
800 (*_field)+=(*f.getField());
804 //----------------------------------------------------------------------
806 Field::operator-= ( const Field& f )
807 //----------------------------------------------------------------------
809 //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
810 //throw CdmathException("Field::operator-= : Field subtraction requires identical meshes");
811 _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
812 if(f.getTypeOfField() != getTypeOfField())
813 throw CdmathException("Field::operator-= : Field subtraction requires identical field types (CELLS, NODES or FACES");
814 if(f.getNumberOfComponents() != getNumberOfComponents())
815 throw CdmathException("Field::operator-= : Field subtraction requires identical number of components");
817 _field->setMesh(f.getField()->getMesh());
818 (*_field)-=(*f.getField());
822 //----------------------------------------------------------------------
824 Field::operator*= ( double s )
825 //----------------------------------------------------------------------
827 int nbComp=getNumberOfComponents();
828 int nbElem=getNumberOfElements();
829 for (int i=0 ; i<nbComp ; i++)
830 for (int j=0 ; j<nbElem; j++)
831 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]*=s;
835 //----------------------------------------------------------------------
837 Field::operator/= ( double s )
838 //----------------------------------------------------------------------
840 int nbComp=getNumberOfComponents();
841 int nbElem=getNumberOfElements();
842 for (int i=0 ; i<nbComp ; i++)
843 for (int j=0 ; j<nbElem; j++)
844 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]/=s;
848 //----------------------------------------------------------------------
850 Field::operator-= ( double s )
851 //----------------------------------------------------------------------
853 int nbComp=getNumberOfComponents();
854 int nbElem=getNumberOfElements();
855 for (int i=0 ; i<nbComp ; i++)
856 for (int j=0 ; j<nbElem; j++)
857 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]-=s;
861 //----------------------------------------------------------------------
863 Field::operator+= ( double s )
864 //----------------------------------------------------------------------
866 int nbComp=getNumberOfComponents();
867 int nbElem=getNumberOfElements();
868 for (int i=0 ; i<nbComp ; i++)
869 for (int j=0 ; j<nbElem; j++)
870 _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]+=s;
874 //----------------------------------------------------------------------
876 Field::writeVTK (std::string fileName, bool fromScratch) const
877 //----------------------------------------------------------------------
879 string fname=fileName+".pvd";
881 double time=_field->getTime(iter,order);
885 ofstream file(fname.c_str()) ;
886 file << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\"><Collection>\n" ;
887 ostringstream numfile;
889 string filetmp=fileName+"_";
890 filetmp=filetmp+numfile.str();
891 string ret=_field->writeVTK(filetmp.c_str()) ;
892 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
893 file << "</Collection></VTKFile>\n" ;
898 ifstream file1(fname.c_str()) ;
900 getline(file1, contenus, '\0');
901 string to_remove="</Collection></VTKFile>";
902 size_t m = contenus.find(to_remove);
903 size_t n = contenus.find_first_of("\n", m + to_remove.length());
904 contenus.erase(m, n - m + 1);
906 ofstream file(fname.c_str()) ;
908 ostringstream numfile;
910 string filetmp=fileName+"_";
911 filetmp=filetmp+numfile.str();
912 string ret=_field->writeVTK(filetmp.c_str()) ;
913 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
914 file << "</Collection></VTKFile>\n" ;
919 //----------------------------------------------------------------------
921 Field::writeCSV ( const std::string fileName ) const
922 //----------------------------------------------------------------------
925 double time=_field->getTime(iter,order);
927 ostringstream numfile;
929 string filetmp=fileName+"_";
930 filetmp=filetmp+numfile.str();
931 filetmp=filetmp+".csv";
932 ofstream file(filetmp.c_str()) ;
933 int dim=_mesh.getSpaceDimension();
935 if (getTypeOfField()==CELLS)
936 nbElements=_mesh.getNumberOfCells();
938 nbElements=_mesh.getNumberOfNodes();
942 int nbCompo=getNumberOfComponents();
944 file << "x " << _field->getName() << endl;
948 for (int i=0;i<nbCompo;i++)
949 file << " " << _field->getName() << " Compo " << i+1 << " "<< getInfoOnComponent(i);
952 for (int i=0;i<nbElements;i++)
954 if (getTypeOfField()==CELLS)
955 file << _mesh.getCell(i).x() ;
957 file << _mesh.getNode(i).x() ;
958 for (int j=0;j<nbCompo;j++)
959 file << " " << getValues()[j+i*nbCompo] ;
964 int nbCompo=getNumberOfComponents();
966 file << "x y " << _field->getName() << endl;
970 for (int i=0;i<nbCompo;i++)
971 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
974 for (int i=0;i<nbElements;i++)
976 if (getTypeOfField()==CELLS)
977 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() ;
979 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() ;
980 for (int j=0;j<nbCompo;j++)
981 file << " " << getValues()[j+i*nbCompo] ;
986 int nbCompo=getNumberOfComponents();
988 file << "x y z " << _field->getName() << endl;
992 for (int i=0;i<nbCompo;i++)
993 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
996 for (int i=0;i<nbElements;i++)
998 if (getTypeOfField()==CELLS)
999 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() << " " << _mesh.getCell(i).z();
1001 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() << " " << _mesh.getNode(i).z();
1002 for (int j=0;j<nbCompo;j++)
1003 file << " " << getValues()[j+i*nbCompo] ;
1010 //----------------------------------------------------------------------
1012 Field::writeMED ( const std::string fileName, bool fromScratch) const
1013 //----------------------------------------------------------------------
1015 string fname=fileName+".med";
1017 MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
1019 MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
1023 operator* (double value , const Field& field )
1025 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1026 int nbComp=field.getNumberOfComponents();
1027 int nbElem=field.getNumberOfElements();
1028 for (int ielem=0 ; ielem<nbElem; ielem++)
1029 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1030 fres(ielem, jcomp)=value*field(ielem, jcomp);
1035 operator* (const Field& field, double value )
1037 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1038 int nbComp=field.getNumberOfComponents();
1039 int nbElem=field.getNumberOfElements();
1040 for (int ielem=0 ; ielem<nbElem; ielem++)
1041 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1042 fres(ielem, jcomp)=value*field(ielem, jcomp);
1046 Field operator/ (const Field& field, double value)
1048 Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1049 int nbComp=field.getNumberOfComponents();
1050 int nbElem=field.getNumberOfElements();
1051 for (int ielem=0 ; ielem<nbElem; ielem++)
1052 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1053 fres(ielem, jcomp)=field(ielem, jcomp)/value;
1058 Field::getValuesOnAllComponents(int elem) const
1060 Vector v(getNumberOfComponents());
1061 for(int i=0;i<getNumberOfComponents();i++)
1062 v[i]=(*this)(elem,i);
1067 Field::getValuesOnComponent(int compo) const
1069 Vector v(getNumberOfElements());
1070 for(int i=0;i<getNumberOfElements();i++)
1071 v[i]=(*this)(i,compo);
1075 std::ostream& operator<<(std::ostream& out, const Field& field )
1077 cout << "Field " << field.getName() << " : " << endl ;
1078 out<< field.getField().retn()->getArray()->repr();