Queries in accessing (I-1)th and (I+1)th (neighboring) Nodes

[kumaresh...@gmail.com]

Hello Foamers, Hope everyone doing good.

Generally, OpenFOAM is designed based on unstructured polyhedral code and the addressing is face-based. My equations are resolved in OpenFOAM by taking into account of differential operators (div,grad,laplacian) in natural way. But, the source term in my custom solver is resolved by node-based approach. Concerning this, I'm hacking finite-volume based code in finite-difference sense as shown below.

//**********************************************************************************//

const Foam::List<int> neigbhorcellI = regionMesh().cellCells(cellI);//gives the neighbours of cell cellI
forAll(neigbhorcellI, I) //Loop inside internal cells
if (regionMesh().cellCentres()[I][0] < regionMesh().cellCentres()[cellI][0])
{
std::cout<<"QB calculation @ cellCentres inside loop"<<std::endl;
QB_[cellI] = K_[I]* ((T_[I] - Tliq) /DX.value()); -----> Equ. (1)
// Equ. (1) solves QB for previous node, lets say (I-1)th node
break ;
}

if (regionMesh().cellCentres()[I][0] > regionMesh().cellCentres()[cellI][0])
{
std::cout<<"QC calculation @ cellCentres inside loop"<<std::endl;
QC_[cellI] = K_[I]* ((Tliq - T_[I]) /DX.value()); -----> Equ. (2)
// I want to solve Equ. (2) for QC at the very next node, lets say (I+1)th node
break ;
}

//**********************************************************************************//

    Am I right in accessing the neighboring nodes --> (I-1)th and (I+1)th nodes as shown above ? By Equ. (1), QB is solved based on (I-1)th node. But QC is not resolved at (I+1)th node (if-loop for QC is not getting executed), here the print statement for QC is not executed in the terminal window. No ERRORS shown while compiling.

Kind of in confusion here, whether I'm right or wrong in accessing neighboring nodes appropriately.

Kindly someone guide me from here please.

Thank you

[penguinitis]

Hello.

The variable I in the forAll loop of neigbhorcellI is the index of neigbhorcellI.

I tested it with code like below and it worked.

        const label cellI = 10; 

        const Foam::List<int> neighbourcellI = mesh.cellCells(cellI);

        forAll(neighbourcellI, I) {

            const label nI = neighbourcellI[I];

            Info << mesh.cellCentres()[nI][0] << " " << mesh.cellCentres()[cellI][0] << endl;

            if (mesh.cellCentres()[nI][0] < mesh.cellCentres()[cellI][0])

            {

                Info<<"QB calculation @ cellCentres inside loop"<<endl;

            }

            if (mesh.cellCentres()[nI][0] > mesh.cellCentres()[cellI][0])

            {

                Info<<"QC calculation @ cellCentres inside loop"<<endl;

            }

        } 

2020年10月9日金曜日 14:35:20 UTC+9 kumaresh...@gmail.com:

[Kumaresh Selvakumar]

Hello penguinitis,

Thank you so much for your prompt response. It really helped me to find a way.

The line below, helped to resolve my problem.

            const label nI = neighbourcellI[I];

I have another query. I would like to know the significance of below line,

        const label cellI = 10;

Here 10 indicates the number of nodes in the mesh ?

Kindly correct me if I am wrong.

Thank you

--
このメールは Google グループのグループ「OpenFOAM」に登録しているユーザーに送られています。
このグループから退会し、グループからのメールの配信を停止するには openfoam+u...@googlegroups.com にメールを送信してください。
このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/ff332845-b661-47f5-9481-80b2fff3c86an%40googlegroups.com にアクセスしてください。

[penguinitis]

> const label cellI = 10;

This is the cell ID.
It is fixed in one cell for testing.

2020年10月12日月曜日 3:33:24 UTC+9 kumaresh...@gmail.com:

[Kumaresh Selvakumar]

Thank you.

> const label cellI = 10;  So 10 is cellId number. It is fixed with one cell for testing as you said.

In my case, the code need to be executed for consecutive cells based on iterations.

//**********************************************************************************//

forAll(T_, cellI)

{
if ((T_[cellI] - Tliq == scalar (0.0)) && (alpha_[cellI] > scalar(0.0)))
{
Info << regionMesh().cellCells(cellI)<< endl; //Print the list of cell index in the neighbourhood

const Foam::List<int> neigbhorcellI = regionMesh().cellCells(cellI);//gives the neighbours of cell cellI
forAll(neigbhorcellI, I) //Loop inside internal cells

{
const label nI = neigbhorcellI[I]; // INDEX necessary to consider the Neigbhour cells                       

Info << "i=" << I << " nI =" << nI << endl;
Info << "C[I]=" << regionMesh().cellCentres()[I] << " C[ni]=" << regionMesh().cellCentres()[nI] << endl;
Info << regionMesh().cellCentres()[nI][0] << " " << regionMesh().cellCentres()[cellI][0] << endl;

if (regionMesh().cellCentres()[nI][0] < regionMesh().cellCentres()[cellI][0])
{
std::cout<<"QB calculation @ cellCentres inside loop"<<std::endl; //1
QB_[cellI] = K_[nI]* ((T_[nI] - Tliq) /DX.value());
//break ;
}

if (regionMesh().cellCentres()[nI][0] > regionMesh().cellCentres()[cellI][0])
{
std::cout<<"QC calculation @ cellCentres inside loop"<<std::endl; //2
QC_[cellI] = K_[nI]* ((Tliq - T_[nI]) /DX.value());
//break ;
}

}

}

//**********************************************************************************//

The above code works well without errors for all consecutive cells. But I didn't achieve my expected results based on this code.

If you found any errors technically on above code, kindly correct me, I would like to learn and improve myself ahead.

Thank you for your precious time.

have a nice day

このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/5c1940cf-c095-40b7-a6a0-d992c6361a9bn%40googlegroups.com にアクセスしてください。

[penguinitis]

>if ((T_[cellI] - Tliq == scalar (0.0)) && (alpha_[cellI] > scalar(0.0)))

Floating point variables rarely have the same value, so I think it's better to have a range for comparison.

For example:

if ((abs(T_[cellI] - Tliq) < 1e-3) && ...

2020年10月12日月曜日 16:41:56 UTC+9 kumaresh...@gmail.com:

[Kumaresh Selvakumar]

Thats valid point.

Thank you. I'll look into it. 

このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/a511197c-8682-453b-ae26-6feabcdb3b6cn%40googlegroups.com にアクセスしてください。

[Kumaresh Selvakumar]

Hello penguinitis,

I have one more query to ask. Previously, the code was implemented to access (I-1) and (I+1) nodes.

//**********************************************************************************//

if (regionMesh().cellCentres()[nI][0] < regionMesh().cellCentres()[cellI][0])    // At (I-1)

{
std::cout<<"QB calculation @ cellCentres inside loop"<<std::endl;
}

if (regionMesh().cellCentres()[nI][0] > regionMesh().cellCentres()[cellI][0])   // At (I+1)

{
std::cout<<"QC calculation @ cellCentres inside loop"<<std::endl;
}

//**********************************************************************************//

T(I) = T(I) - X *[ (T(I-1) - T(I+1) ] / 2;

As per this T(I) equation, how to access both (I-1) and (I+1) by one equation ?

Any ideas ?

Thank you

Regards

Kumaresh

[Kumaresh Selvakumar]

Hello penguinitis,

For the question above, will this answer below be an appropriate one ?

//**********************************************************************************//

if (regionMesh().cellCentres()[nI][0] < regionMesh().cellCentres()[cellI][0])    // At (I-1)

{
T_[nI] = dummy1_[cellI];
}

if (regionMesh().cellCentres()[nI][0] > regionMesh().cellCentres()[cellI][0])   // At (I+1)

{
T_[nI] = dummy2_[cellI];
}

//**********************************************************************************//

T(I) = T(I) - X *[ (T(I-1) - T(I+1) ] / 2;

==> T(I) = T(I) - X *[ (dummy1_[cellI] - dummy2_[cellI]) ] / 2;

Thank you

Regards

Kumaresh

[penguinitis]

Do you mean that you consider the center of the cell to be a node?

If you want to calculate the value of a cell from the left and right cells, you can save the indexes of the left and right cells.

        const Foam::List<int> neighbourcellI = mesh.cellCells(cellI);

        label nI0 = -1, nI1 = -1;

        forAll(neighbourcellI, I) {
            const label nI = neighbourcellI[I];

            if (mesh.cellCentres()[nI][0] < mesh.cellCentres()[cellI][0])
            {

                nI0 = nI;

            }
            if (mesh.cellCentres()[nI][0] > mesh.cellCentres()[cellI][0])
            {

                nI1 = nI;
            }
            if (nI0 > 0 && nI1 > 0 ) break;
        }
        Info << nI0 << " " << nI1 << endl;

T(I) = T(I) - X *[ (T(I-1) - T(I+1) ] / 2;

==> T[cellI] = T[cellI] - X *(T[nI0] - T[nI1] )/ 2;

2020年10月14日水曜日 23:23:40 UTC+9 kumaresh...@gmail.com:

[Kumaresh Selvakumar]

Hello penguinitis,

Thank you so much for your valuable feedback once again.

Yes, you are right, I consider the centre of the cell.

T[cellI] = T[cellI] - X *(T[nI0] - T[nI1] )/ 2;

Here this T loop is executed well without any compilation errors.

However, I have some queries to ask:

1) label nI0 = -1, nI1 = -1; --> Here, may I know why the indexes are declared as -1.

2) if (nI0 > 0 && nI1 > 0 ) break; --> What is the significance of this loop ?

Thank you

Regards

Kumaresh

このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/44b82d26-a4e3-455d-80b3-35fd2d72ebffn%40googlegroups.com にアクセスしてください。

[penguinitis]

1) 2) These are for exiting the loop if both adjacent cells are found.

2) was a little wrong. The following is correct.

if (nI0 >= 0 && nI1 >= 0 ) break;

If your model is one-dimensional, you don't need these.

2020年10月16日金曜日 20:07:55 UTC+9 kumaresh...@gmail.com:

[Kumaresh Selvakumar]

Hello penguinitis,

         Thank you for your email and time.

Your thoughts helped me a lot.

Thank you

Regards

Kumaresh

このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/3e7a3636-2490-4dd5-afed-14df8eaae893n%40googlegroups.com にアクセスしてください。

[kumaresh...@gmail.com]

Hello penguinitis,

I have one more query to ask.

//**********************************************************************************//

forAll(T_, cellI)   
{

    if (T_[cellI] = Tliq)   // Access at (I)th cell
    {

          forAll(neigbhorcellI, I)
          {
          const label nI = neigbhorcellI[I];                                            
          if (regionMesh().cellCentres()[nI][0] < regionMesh().cellCentres()[cellI][0]) // This should access (I-1) cell
          {   

                T_[nI] = T_[nI] - XX *(T_[I-2] - T[I]) / 2;  // Calculation at (I-1) is based on (I-2) and (I) --> Eq. (1)

          }

   }

} 

//**********************************************************************************//

Here,

T(I) = Tliq (fixed value, no issues here)

T(I-1) = T_[nI] (need to update at (I-1) node based on if loop using Eq. (1), no issues here)

T(I-2) = How to calculate the value at (I-2) node ?

I'm hereby to ask you query regarding, as how to calculate the value at (I-2) node ?

Kindly share your ideas. Thank you.

Regards

Kumaresh

[penguinitis]

Hello.

How about creating tables of the index for the neighbour cell in advance?

Assuming a 1D mesh:

//////////

Foam::Map<label> niTable0, niTable1;  // index tables

forAll(T_, cellI) {

        const Foam::List<label> neighbourcellI = mesh.cellCells(cellI);

        label nI0 = -1, nI1 = -1;                                                
                                                                                 
        forAll(neighbourcellI, I) {                                              
            const label nI = neighbourcellI[I];                                  
                                                                                 

            Info << mesh.cellCentres()[nI][0] << " " << mesh.cellCentres()[cellI][0] << endl;
                                                                                 
            if (mesh.cellCentres()[nI][0] < mesh.cellCentres()[cellI][0])        
            {

                nI0 = nI;                                                        
            }                                                                    
                                                                                 
            if (mesh.cellCentres()[nI][0] > mesh.cellCentres()[cellI][0])        
            {
                nI1 = nI;                                                        
            }                           
        }                                                                        
                                                                                 

        niTable0.insert(cellI, nI0);                                             
        niTable1.insert(cellI, nI1);                                             
                                                                                 
        Info << cellI << " " << niTable0[cellI] << " " << niTable1[cellI] << endl;                                   
        Info << "check: " << cellI << " " << nI0 << " " << nI1 << endl;

}

forAll(T_, cellI)   
{

    if (T_[cellI] = Tliq)

    {

        ...

//////////

nI = niTable0[cellI]

nI2 = niTable0[nI]

T(I-1) = T_[nI]

T(I-2) = T_[nI2]

2020年11月4日水曜日 19:20:22 UTC+9 kumaresh...@gmail.com:

[kumaresh...@gmail.com]

Thank you penguinitis for your clarified thoughts again.

The above technique about creating a table of cell indexes is a good hit for my work. Thank you again ^^

I got some errors, when I tried implementing the following lines:

nI = niTable0[cellI]

nI2 = niTable0[nI]

ERRORS:

 >> ‘nI’ was not declared in this scope

 >> ‘nI2’ was not declared in this scope

The thoughts are well clear, that 'nI2' can be applied for T(I-2). But somewhere the declaration of 'nI2' is a issue.

Kindly check it when you are free please.

Thank you

Have a nice day ^^

[kumaresh...@gmail.com]

Hello penguinitis,

     I have tried as below,

 scalar kum = niTable0[cellI];
 scalar nI2 = niTable0[kum];

 

I declared as scalar variables and tried fixing it appropriately.

T_[nI1] = T_[nI1] - XX *(T_[nI2] - Tliq) / 2;

No compilation error here. Is this okay ?

Kindly share your thoughts.

Thank you

[penguinitis]

Hello Kumaresh,

nI = niTable0[cellI]
nI2 = niTable0[nI]

This part was meant to be an explanation, not a code.

These types are "label".

label nI = niTable0[cellI];
label nI2 = niTable0[nI];

2020年11月7日土曜日 18:08:27 UTC+9 kumaresh...@gmail.com:

[kumaresh...@gmail.com]

Thank you for your thoughts again ^^

[Jaime Andres Castañeda Barbosa]

Good morning dears,

I am trying to run a case but in order to run I need to have installed swak4Foam in OpenFOAM. 

I want to know if anyone has any videos on how to install swak4Foam (GroovyBC) in OpenFoam ?. 

I'm not able to install on my Ubuntu virtual machine (Ubuntu 7.5.0-3ubuntu1~18.04) 7.5.0. the version is swak4Foam 0.3.2.

https://openfoamwiki.net/index.php/Installation/swak4Foam/Installing_On/Ubuntu

 

Thank you for your time and consideration,

[kumaresh...@gmail.com]

Hello penguinitis,

Hope you are doing good. Long before, I got the help from you for accessing (I), (I-1) and (I-2) nodes. I have some query here to follow.

Let me explain the code first.

Assuming a 1D mesh:

Foam::Map<label> niTable0, niTable1;

forAll(T_, cellI)

{

/////////////////

        const Foam::List<label> neighbourcellI = mesh.cellCells(cellI);

        label nI0 = -1, nI1 = -1;                                                
                                                                                 
        forAll(neighbourcellI, I) {                                              
            const label nI = neighbourcellI[I];                                 

            if (mesh.cellCentres()[nI][0] < mesh.cellCentres()[cellI][0])        
            {              

                 nI0 = nI; 
            }                                                                    
                                                                                 
            if (mesh.cellCentres()[nI][0] > mesh.cellCentres()[cellI][0])        
            {
                nI1 = nI;                                                        
            }                           
        }                                                                                                                                                 

        niTable0.insert(cellI, nI0);                                             
        niTable1.insert(cellI, nI1);                                            

/////////////

            label nI = niTable0[cellI];    //T(I-1)
            label nI2 = niTable0[nI];    //T(I-2)

        if ((abs(T_[cellI] - Tliq) < 1e-2)))
       {

            T_[nI0] = T_[nI0] - XX *(T_[nI2] - Tliq) / 2;        ----> EQU. (1)

            //T(I)=T(I)- XX*(T(I-1)-T(I+1))/2;   

        }

}

While running the case file, the following ERROR occurs at EQU. (1):

//*************************************************************//
--> FOAM FATAL ERROR:
index -1 out of range 0 ... 36

    From function UList<T>::checkIndex(const label)
    in file /home/kummi/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/UListI.H at line 109.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-2.1.1/src/OSspecific/POSIX/printStack.C:201
#1  Foam::error::abort() at ~/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/error.C:249
#2  Foam::Ostream& Foam::operator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) at ~/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/errorManip.H:85 (discriminator 3)
#3  Foam::UList<double>::checkIndex(int) const at ~/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/UListI.H:113
#4  Foam::UList<double>::operator[](int) at ~/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/UListI.H:168

//*************************************************************//

The above error says that @UlistI.H at line 109,        

FatalErrorIn("UList<T>::checkIndex(const label)")
            << "index " << i << " out of range 0 ... " << size_-1
            << abort(FatalError);

So, what I feel the reason of an error, could be due to nI2 (I-2)th node. May be because it is not initialized as how nI0 and nI1 set to -1 //label nI0 = -1, nI1 = -1;

I'm not sure about my guess.

Kindly help me out here to resolve the error and correct me if I'm wrong.

Thank you ^^

On Monday, November 9, 2020 at 3:15:24 PM UTC+5:30 penguinitis wrote: