# Integration(math) issue

 Register Blogs Members List Search Today's Posts Mark Forums Read

November 9, 2015, 16:50
Integration(math) issue
#1
Member

Sandra
Join Date: Oct 2014
Posts: 58
Rep Power: 4
Hello Foamers!
This code might look familiar but I have added a few lines that require an integration that isn't going so well. Anyhoo, here's the code

Code:
forAll(x3, celli)
{
for (count=0; count<iter; count++)
{
x3 = (x1 + x2)/2;

jc = i*(Foam::exp(8.1*x))-(Foam::exp(-2.7*x))

f1 = (Foam::exp(8.1*x1))-(Foam::exp(-2.7*x1))-(jc/i);

f3 = (Foam::exp(8.1*x3))-(Foam::exp(-2.7*x3))-(jc/i);

if(neg(f1[celli] *f3[celli]))
{
x2 = x3;
}
else if((mag(x1[celli] - x2[celli]) < 1e-7)||(f3[celli]==0)||(jc[celli].integrate(0,3.6e9)==2e4))
{
root =x3;
break;
}
else
{
x1= x3;
}
}

Info<<"Root is"<<endl;
Info<<root<<endl;
}
jc is declared as a volScalarField in my createFields and here's the error that I get (I understand that I do not have a double ie I understand what the error means but I do not know how to proceed.

Quote:
 error: request for member ‘integrate’ in ‘jc.Foam::GeometricField::.Foam:imensionedField ::.Foam::Field::< anonymous>.Foam::List::.Foam::U List:perator[](celli)’, which is of non-class type ‘double’
Sandra

 November 10, 2015, 04:20 #2 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,691 Rep Power: 28 Hi, 1. Code: f3[celli]==0 in general it is not a good idea to compare doubles with ==. There are two functions: Code: inline bool equal(const Scalar& s1, const Scalar& s2) { return mag(s1 - s2) <= ScalarVSMALL; } inline bool notEqual(const Scalar s1, const Scalar s2) { return mag(s1 - s2) > ScalarVSMALL; } 2. if jc is volScalarField then jc[celli] is just scalar (i.e. simple number), so And again this Code: jc[celli].integrate(0,3.6e9)==2e4 will be FALSE 99,9% of the time.

November 10, 2015, 04:51
#3
Member

Sandra
Join Date: Oct 2014
Posts: 58
Rep Power: 4
Quote:
 Originally Posted by alexeym 2. if jc is volScalarField then jc[celli] is just scalar (i.e. simple number), so And again this Code: jc[celli].integrate(0,3.6e9)==2e4 will be FALSE 99,9% of the time.
Would that formula work even if C is not a linear function per say (and I don't know its form either)? There's no method in OpenFOAM to get the area under a curve during simulation for volScalarFields like there's for doubles ie the integrate that I tried to use?

Lastly, out of curiosity, why would
Code:
jc[celli].integrate(0,3.6e9)==2e4
be false 99.9% of the time?

November 10, 2015, 05:10
#4
Senior Member

Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,691
Rep Power: 28
If, as you said, jc is volScalarField, then jc[celli] can not be "...not a linear function per say...". It is a number and on a given time step it is in general even constant. I guess, you need to explain what you are trying to achieve with jc[celli].integrate(...). Right now what you are trying to do is like this:

Code:
(3.1415).integrate(0, 3e9)
Why do you need to integrate constant? What is integration variable?

Quote:
 Lastly, out of curiosity, why would ... be false 99.9% of the time?
https://isocpp.org/wiki/faq/newbie#floating-point-arith
https://en.wikipedia.org/wiki/Floati...uracy_problems
https://duckduckgo.com/?q=floating+point+comparison

 November 10, 2015, 05:27 #5 Member   Sandra Join Date: Oct 2014 Posts: 58 Rep Power: 4 sorry, my bad.... jc in paraview (at each time step) is a curve like a straight up witch hat and I want to iterate until the area under that curve equals a constant that I give it (in this case, 20,000). So, I agree, jc[celli] wouldn't work here but I wasn't thinking of that, unfortunately, when I tried it after "jc.integrate" gave me an error along the lines of "member doesn't have function integrate"

 November 10, 2015, 05:37 #6 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,691 Rep Power: 28 Well, let us start from the beginning... What is jc? In the first post you wrote, that it is volScalarField, i.e. something that varies in VOLUME, then you have written, that it is "...a curve like a straight up witch hat...", i.e. someting that varies just in ONE direction (surely there are 3D curves, yet they usually can be reduced to function of one variable). After you define jc, I think, it is easier to invent a method for calculation your "integral".

 November 10, 2015, 05:41 #7 Member   Sandra Join Date: Oct 2014 Posts: 58 Rep Power: 4 it is defined as a volScalarField and calculated from the equation in the initial post BUT my case is 2D ie I have front and back patches empty so technically it's an areaScalarField if that were a thing in OpenFOAM and that's why I talked about a 2d curve. Does that help? Sandra

 November 10, 2015, 05:58 #8 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,691 Rep Power: 28 OK, your case is 2D, so is jc a surface? So you have and would like to calculate where xa, xb, ya, yb are the corners of your simulation domain (so it would be volume under the curve and in simplest form it is a sum of jc[celli]*mesh.V()[celli]/)? Or you would like to take values of jc along the line L (technically reducing surface to curve) and them integrate 1D function along L?

November 10, 2015, 06:03
#9
Member

Sandra
Join Date: Oct 2014
Posts: 58
Rep Power: 4
Quote:
 Originally Posted by alexeym OK, your case is 2D, so is jc a surface? So you have and would like to calculate where xa, xb, ya, yb are the corners of your simulation domain (so it would be volume under the curve and in simplest form it is a sum of jc[celli]*mesh.V()[celli]/)? Or you would like to take values of jc along the line L (technically reducing surface to curve) and them integrate 1D function along L?

Many thanks, BTW

 November 10, 2015, 06:30 #10 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,691 Rep Power: 28 There are several possibilities: 1. Your mesh is static regular orthogonal, your line is parallel to axis. You construct cell set (cells, where you line passes) yourself in the beginning of the simulation, starting from xa (using, for example, findCell method of mesh) and then iterating through the neighbors and selecting cells in required direction (or if mesh is not large, you can use findCell for every point of your line, as you are doing it just once). After in simplest case you iterate through the labels, take value of jc in the cell, take width of the cell and multiply. 2. You mesh is static irregular, your line is defined as a*x + b. You can use polyLineSet set class to construct cell set. 3. You mesh moves, is irregular, your line is defined as a*x + b. You reconstruct cell set using polyLineSet after every mesh topological change.

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post scnsc OpenFOAM Installation 3 August 20, 2014 07:59 Philip_C ANSYS 11 August 29, 2013 06:44 prasant OpenFOAM Running, Solving & CFD 17 March 16, 2013 03:00 Vijay FLUENT 0 April 6, 2012 13:35 appu FloEFD, FloWorks & FloTHERM 1 May 22, 2011 08:27

All times are GMT -4. The time now is 07:46.