[deal.II] Re: Regarding configuring with UMFPACK for Quasi-Static Strain Compressible Elasticity program

2020-09-12 Thread Simon Sticko
Hi,
You'll need to call cmake with the flag

-DDEAL_II_WITH_UMFPACK = ON

> Also, do I need to reinstall dealii everytime I realise that I have not 
configured dealii with a particular package?

Yes. However, if you need several of the optional dependencies, there is a 
script that can make the installation process easier here:
https://github.com/dealii/candi

Best,
Simon

On Saturday, September 12, 2020 at 1:42:29 PM UTC+2, Animesh Rastogi IIT 
Gandhinagar wrote:
>
> Attached is the detailed.log file. 
> On Saturday, September 12, 2020 at 4:39:21 PM UTC+5:30 Animesh Rastogi IIT 
> Gandhinagar wrote:
>
>> Hello dealii Community, 
>>
>> I wish to run the program - Quasi-Static Finite-Strain Compressible 
>> Elasticity. 
>>
>> 
>>
>> This requires the dealii to be configured with UMFPACK. For this, I 
>> installed the SuiteSparse with UMFPACK, compiled it, and then reconfigured 
>> and reinstalled dealii using the command -DUMFPACK_DIR=/path/to/umfpack 
>> (changing the path to the compiled version of UMFPACK in my PC) with cmake.
>>
>> However, I am still getting the following error when I am running the 
>> Finite Strain program.
>>
>> An error occurred in line <635> of file 
>>  
>> in function
>> void dealii::SparseDirectUMFPACK::factorize(const Matrix&) [with 
>> Matrix = dealii::SparseMatrix]
>> The violated condition was: 
>> false
>> Additional information: 
>> To call this function you need UMFPACK, but you configured deal.II 
>> without passing the necessary switch to 'cmake'. Please consult the 
>> installation instructions in doc/readme.html.
>>
>> Could someone please let me know how should I sort it out?
>>
>> Also, do I need to reinstall dealii everytime I realise that I have not 
>> configured dealii with a particular package?
>>
>> Thanks!
>>
>> Animesh
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/61ee8470-ed18-42fa-ade4-b74d208a0809o%40googlegroups.com.


[deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder

2020-07-01 Thread Simon Sticko
Hi,

> and the code seems to be working fine! Is there anything that I am 
> missing in your opinion?

No. If you're happy, I'm happy.

Best,
Simon

On Wednesday, July 1, 2020 at 7:51:11 AM UTC+2, Samuel Rodriguez wrote:
>
> Hi Simon,
>
> I am a slow learner but you were correct. Here is the first image of the 
> function plotted along the radius. As you can see there is only noise but 
> value of the solution is constant. I can also change the value of the 
> solution by changing the value of the boundary condition as you can also 
> see below:
>
> [image: Screen Shot 2020-06-30 at 2.20.11 PM.png]
>
> [image: Screen Shot 2020-06-30 at 10.22.07 PM.png]
>
> I was able to change my boundary condition to something ridiculous like:
>
> class NeumannBoundaryCondition : public Function<3>
> {
> public:
>   double
>   value(const Point<3> ,  const unsigned int component = 0) const
>   {
> double test = 1;
> for(int i = 0; i < 3; i++)
>   test += 20 + std::pow(point(i), 20) + std::pow(point(i), 10) + 
> exp(point(i));
> std::cout << test << std::endl;
> return test;
>   }
> };
>
>
>
> The reason this was not initially working was because when I implemented 
> the boundary condition below:
>
>   for(unsigned int face_number = 0; face_number < 
> GeometryInfo<3>::faces_per_cell; ++face_number)
>   {
>if(cell->face(face_number)->at_boundary() && 
> (cell->face(face_number)->boundary_id() == 0))
> {
>  fe_face_values.reinit(cell, face_number);
>  for (unsigned int q_point = 0; q_point < n_face_q_points; 
> ++q_point)
>   {
>double neumannvalue = 
> boundary_condition.value(fe_face_values.quadrature_point(q_point));
>for (unsigned int i = 0; i < dofs_per_cell; ++i)
> {
>  const unsigned int component_i = 
> fe.system_to_component_index(i).first;
>  cell_rhs(i) += (neumannvalue *
>  fe_face_values.shape_value(i, q_point) *
>  fe_face_values.JxW(q_point));
> }
>   }
> 
>}
>  }
>
>
> I never added it the right hand side (I think that I am using that 
> terminology correctly). I added the previous lines of code right above these
>
>   cell->get_dof_indices(local_dof_indices);
>
>   for (unsigned int i = 0; i < dofs_per_cell; ++i)
> for (unsigned int j = 0; j < dofs_per_cell; ++j)
>   system_matrix.add(local_dof_indices[i],
> local_dof_indices[j],
> cell_matrix(i, j));
>   for (unsigned int i = 0; i < dofs_per_cell; ++i)
> system_rhs(local_dof_indices[i]) += cell_rhs(i);
>
>
> and the code seems to be working fine! Is there anything that I am missing 
> in your opinion?
>
> Best,
>
> Samuel Rodriguez
>
>
>
>
>
>
>
>
>
>
>
> On Tuesday, June 30, 2020 at 3:33:25 AM UTC-7, Simon Sticko wrote:
>>
>> Hi,
>> From the figure you can't really tell whether the boundary condition is 
>> fulfilled or not. If you want to check if the normal derivative is 
>> reasonable (almost equal to 0) at the boundary, the easiest is probably to 
>> plot the solution on a line in the radial direction. This can be done by 
>> using
>> Filters -> Data Analysis -> Plot Over Line in Paraview.
>>
>> I don't quite understand what you mean when you say that you have tried 
>> changing the value of component. Component is intended to be used for 
>> vector valued problems. If you are solving Poisson, it doesn't make sense 
>> to set it to something different than 0. 
>>
>> Note that, right now you are calling the function 
>> NeumannBoundaryCondition::value with a single argument (of type Point<3>). 
>> This means that you get the default argument for component: component=0.
>>
>> Best,
>> Simon
>>
>> On Tuesday, June 30, 2020 at 2:59:20 AM UTC+2, Samuel Rodriguez wrote:
>>>
>>> Hi Simon,
>>>
>>> It seems I might not have understood my own question. You answered the 
>>> question I was trying to ask. Would you be able to help me out with one 
>>> more thing? I followed your advice but the boundary conditions do not seem 
>>> bedo not seem to be getting the correct boundary conditions, so I don't 
>>> think they are being applied correctly. Here is the simple piece of code 
>>> that I wrote:
>>>
>>> class NeumannBoundaryCondition : public Function<

[deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder

2020-06-30 Thread Simon Sticko
Hi,
>From the figure you can't really tell whether the boundary condition is 
fulfilled or not. If you want to check if the normal derivative is 
reasonable (almost equal to 0) at the boundary, the easiest is probably to 
plot the solution on a line in the radial direction. This can be done by 
using
Filters -> Data Analysis -> Plot Over Line in Paraview.

I don't quite understand what you mean when you say that you have tried 
changing the value of component. Component is intended to be used for 
vector valued problems. If you are solving Poisson, it doesn't make sense 
to set it to something different than 0. 

Note that, right now you are calling the function 
NeumannBoundaryCondition::value with a single argument (of type Point<3>). 
This means that you get the default argument for component: component=0.

Best,
Simon

On Tuesday, June 30, 2020 at 2:59:20 AM UTC+2, Samuel Rodriguez wrote:
>
> Hi Simon,
>
> It seems I might not have understood my own question. You answered the 
> question I was trying to ask. Would you be able to help me out with one 
> more thing? I followed your advice but the boundary conditions do not seem 
> bedo not seem to be getting the correct boundary conditions, so I don't 
> think they are being applied correctly. Here is the simple piece of code 
> that I wrote:
>
> class NeumannBoundaryCondition : public Function<3>
> {
> public:
>   double
>   value(const Point<3> ,  const unsigned int component = 0) const
>   {
> return component;
>   }
> };
>
>
> This neumann boundary condition is only meant to apply a constant value of 
> zero on the sides. Here is how I implemented it:
>
> QGauss<3> quadrature_formula(fe.degree + 1);
> QGauss<3 - 1> face_quadrature_formula(fe.degree + 1);
> FEFaceValues<3> fe_face_values(fe,
>face_quadrature_formula,
>update_values | update_quadrature_points |
>update_normal_vectors |
>update_JxW_values);const NeumannBoundaryCondition 
> boundary_condition;
>
> const unsigned int n_face_q_points = face_quadrature_formula.size();
>
>
> for(unsigned int face_number = 0; face_number < GeometryInfo<3>::
> faces_per_cell; ++face_number)
> {
>   if(cell->face(face_number)->at_boundary() && (cell->face(face_number)->
> boundary_id() == 0))
>   {
>fe_face_values.reinit(cell, face_number);
>for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) {
>  const double neumannvalue = boundary_condition.value(fe_face_values.
> quadrature_point(q_point));
>  for (unsigned int i = 0; i < dofs_per_cell; ++i) {
>const unsigned int component_i = fe.system_to_component_index(i).
> first;
>cell_rhs(i) += (neumannvalue *
>fe_face_values.shape_value(i, q_point) *
>fe_face_values.JxW(q_point));
>}
>  }
>  
>   }
> }
>
>
> This seems correct to me, but the code is showing a cyldiner that looks 
> like this:
>
>
> [image: Screen Shot 2020-06-29 at 5.50.19 PM.png]
> which doesn't appear to be the right solution. Even when I change the 
> value of component to something ridiculous such as 10,000,000 the solution 
> on the cylinder does not change. Would you be able to tell me what I am 
> doing wrong in the implementation?
>
> - Samuel
>
>
>
> On Friday, June 26, 2020 at 2:21:19 AM UTC-7, Simon Sticko wrote:
>>
>> Hi,
>> I'm not sure I understand your question correctly, but if you want to 
>> impose an inhomogeneous Neumann boundary condition you would typically 
>> write a small class representing the boundary condition you want:
>>
>> class NeumannBoundaryCondition : public Function<3>
>> {
>> public:
>>   double
>>   value(const Point<3> , const unsigned int component = 0) const
>>   {
>> // Implement the boundary condition you want here.
>> return 1;
>>   }
>> };
>>
>> Then just use this when you assemble:
>>
>> const NeumannBoundaryCondition boundary_condition;
>>
>> const Point<3>  = fe_face_values.quadrature_point(q_point);
>>
>> const double neumann_value = boundary_condition.value(point);
>>
>> cell_rhs(i) += (neumann_value *
>> fe_face_values.shape_value(i, q_point) *
>> fe_face_values.JxW(q_point));
>>
>>
>> Does that answer your question?
>>
>> Best,
>> Simon
>>
>> On Friday, June 26, 2020 at 6:02:48 AM UTC+2, Samuel Rodriguez wrote:
>>>
>>> Hello Everyone,
>>>
>>> I am a masters student barely getting into working with the very heavy 
>>> theoreti

[deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder

2020-06-26 Thread Simon Sticko
Hi,
I'm not sure I understand your question correctly, but if you want to 
impose an inhomogeneous Neumann boundary condition you would typically 
write a small class representing the boundary condition you want:

class NeumannBoundaryCondition : public Function<3>
{
public:
  double
  value(const Point<3> , const unsigned int component = 0) const
  {
// Implement the boundary condition you want here.
return 1;
  }
};

Then just use this when you assemble:

const NeumannBoundaryCondition boundary_condition;

const Point<3>  = fe_face_values.quadrature_point(q_point);

const double neumann_value = boundary_condition.value(point);

cell_rhs(i) += (neumann_value *
fe_face_values.shape_value(i, q_point) *
fe_face_values.JxW(q_point));


Does that answer your question?

Best,
Simon

On Friday, June 26, 2020 at 6:02:48 AM UTC+2, Samuel Rodriguez wrote:
>
> Hello Everyone,
>
> I am a masters student barely getting into working with the very heavy 
> theoretical foundations of using deal.ii. I have been working on a research 
> project that uses a code that creates a cylinder with Dirichlet boundary 
> conditions on the hull/bottom faces and a Neumann boundary condition on the 
> top face Recently I have been trying to change the Dirichlet boundary 
> condition on the hull of the cylinder to a Neumann boundary condition. Our 
> code can be found here inside of QuasistaticBrownianThermalNoise.cpp:
>
> Numerical Coating Thermal Noise 
> 
>
> I have applied the boundary condition but it does not seem to be working. 
> So what I did instead was go to the examples in the tutorial section to 
> learn ow to apply non-homogenous boundary conditions on a cylinder. I have 
> solved the laplacian on the this cylinder using step 3. I have tried to use 
> step 7 to implement non-homogeneous boundary conditions but step 7 requires 
> using a known solution. I would like to be able to set the derivative of 
> the function on the hull to zero and the derivative on the top face to some 
> known function. A Dirichlet boundary condition would be set on the bottom 
> face. The assemble function is down below. My question is how exactly do I 
> apply the Neumann boundary condition to the faces of the cylinder. In the 
> code below I have comments that explain my question in more detail. My 
> question begins where I loop over the faces of the cylinder. However, I 
> included the entire assembly function for completeness in case it was 
> needed. I also attached the whole cpp file in case that was also needed. 
> Any help would be greatly appreciated. If I can provide any more 
> information I would be happy to do so.
>
> QGauss<3> quadrature_formula(fe.degree + 1);
> QGauss<3 - 1> face_quadrature_formula(fe.degree + 1);
>
> FEValues<3> fe_values(fe,
> quadrature_formula,
> update_values | update_gradients | update_JxW_values);
> FEFaceValues<3> fe_face_values(fe,
> face_quadrature_formula,
> update_values | update_quadrature_points |
> update_normal_vectors |
> update_JxW_values);
>
> const unsigned int dofs_per_cell = fe.dofs_per_cell;
> const unsigned int n_q_points = quadrature_formula.size();
> const unsigned int n_face_q_points = face_quadrature_formula.size();
>
> FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
> Vector cell_rhs(dofs_per_cell);
>
> std::vector local_dof_indices(dofs_per_cell);
>
> for (const auto  : dof_handler.active_cell_iterators())
> {
> fe_values.reinit(cell);
>
> cell_matrix = 0;
> cell_rhs = 0;
>
> for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
> {
> for (unsigned int i = 0; i < dofs_per_cell; ++i)
> for (unsigned int j = 0; j < dofs_per_cell; ++j)
> cell_matrix(i, j) +=
> (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
> fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
> fe_values.JxW(q_index)); // dx
> for (unsigned int i = 0; i < dofs_per_cell; ++i)
> cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
> 1 * // f(x_q)
> fe_values.JxW(q_index)); // dx
> }
> cell->get_dof_indices(local_dof_indices);
>
> for (unsigned int i = 0; i < dofs_per_cell; ++i)
> for (unsigned int j = 0; j < dofs_per_cell; ++j)
> system_matrix.add(local_dof_indices[i],
> local_dof_indices[j],
> cell_matrix(i, j));
> for (unsigned int i = 0; i < dofs_per_cell; ++i)
> system_rhs(local_dof_indices[i]) += cell_rhs(i);
>
> // I know that the hull of the cylinder has boundary id 0, the top face 
> has boundary id 1, and the bottom of the face has boundary id 2.
> for(unsigned int face_number = 0; face_number < 
> GeometryInfo<3>::faces_per_cell; 
> ++face_number)
> {
> if(cell->face(face_number)->at_boundary() && 
> (cell->face(face_number)->boundary_id() == 0))
> {
> fe_face_values.reinit(cell, face_number);
> // If we come in here we have found a face that belongs to the boundary 
> condtion of the hull
> // I know that I am supposed to do something like the code in green below, 

Re: [deal.II] How to transform points on faces from unit cells to real cells?

2020-06-22 Thread Simon Sticko
I just want to clarify that Dougs suggestion also works,
if we use the QProjector class to lift the points
on the face from R^(dim-1) to R^dim before we transform them:

// Create a dummy quadrature from support points on face.
const Quadrature dummy_face_quadrature(
fe.get_unit_face_support_points());

// Lift support points on face from R^(dim-1) to R^dim
const Quadrature dummy_quadrature =
QProjector::project_to_face(dummy_face_quadrature, face_no);

// Now we can use transform_unit_to_real_cell
const Point real_point = mapping.transform_unit_to_real_cell(
cell, dummy_quadrature.point(quadrature_points_no));

Depending on how your code looks, this might be an easier solution.

Best,
Simon


On Monday, June 22, 2020 at 11:17:35 PM UTC+2, Lex Lee wrote:
>
> Doug, thanks for sharing me your ideas and the discussion. -Lex
>
>
> On Monday, June 22, 2020 at 2:13:25 PM UTC-7, Doug Shi-Dong wrote:
>>
>> Agreed. Sorry for the confusion. Looked back at my code and it's probably 
>> the most straightforward way.
>>
>> For some other unrelated reason I didn't want to use FEValues and 
>> directly deal with the Mapping class.
>>
>> Doug
>>
>> On Monday, June 22, 2020 at 5:05:57 PM UTC-4, Simon Sticko wrote:
>>>
>>> Hi,
>>>
>>> Lex, given your previous question:
>>>
>>> https://groups.google.com/forum/#!topic/dealii/xghVE7Km78o
>>>
>>> If you are already using an FEFaceValues object with a dummy-quadrature 
>>> (created from FiniteElement::get_unit_face_support_points())
>>> then you can use one of the following functions on FEFaceValues to get 
>>> the location of the support points in real space:
>>>
>>> quadrature_point()
>>> get_quadrature_points()
>>>
>>> Best,
>>> Simon
>>>
>>> On Monday, June 22, 2020 at 10:53:10 PM UTC+2, Lex Lee wrote:
>>>>
>>>> Hello Doug,
>>>>
>>>>
>>>> Thanks for your help.
>>>>
>>>> However, I need to say, I just had done exactly what you described 
>>>> before I posted this question.
>>>>
>>>> "FiniteElement::get_unit_face_support_points()” returns Points;
>>>>
>>>> “transform_unit_to_real_cell()” requires Point as input.
>>>>
>>>> The dimensions of points in cells is different from that on faces. 
>>>> That’s one of the reasons why I posted this question.
>>>>
>>>>
>>>> If I misunderstand you, please kindly let me know.
>>>>
>>>>
>>>>
>>>> Regards,
>>>> Lex 
>>>>
>>>>
>>>> On Jun 21, 2020, at 7:05 PM, Doug Shi-Dong  wrote:
>>>>
>>>> Hello Lex,
>>>>
>>>> transform_unit_to_real_cell() takes a point as unit cell's point an 
>>>> input and returns the physical point location. Therefore, you can give the 
>>>> unit face support points from the 
>>>> FiniteElement::get_unit_face_support_points() (assuming your FiniteElement 
>>>> has support points). And pass it to transform_unit_to_real_cell().
>>>>
>>>> Doug
>>>>
>>>> On Sunday, June 21, 2020 at 8:10:30 PM UTC-4, Lex Lee wrote:
>>>>>
>>>>> Hello Doug,
>>>>>
>>>>> Thanks for your kind help.
>>>>>
>>>>> I am trying to understand your suggestion in the original email.
>>>>>
>>>>> Are you suggesting me using the member function 
>>>>> "transform_unit_to_real_cell()" to get all support points in a unit cell? 
>>>>> In the step, sort out the face points I need? Am I correct?
>>>>>
>>>>>
>>>>> Regards,
>>>>>
>>>>> Lex
>>>>>
>>>>> On Sunday, June 21, 2020 at 5:00:35 PM UTC-7, Doug Shi-Dong wrote:
>>>>>>
>>>>>> Hello Lex,
>>>>>>
>>>>>> You were close.
>>>>>>
>>>>>>
>>>>>> https://www.dealii.org/current/doxygen/deal.II/classMapping.html#ae5df63553eb8ed170c3b90524853dd48
>>>>>>
>>>>>> The project_real_point_to_unit_point_on_face() is useful is you take 
>>>>>> a "volume" point within the cell. Since you are mapping from unit to 
>>>>>> real, 
>>>>>> you would always know whether your point is on the face or not.
>>>>>>
>>>>>> Therefore, tr

Re: [deal.II] How to transform points on faces from unit cells to real cells?

2020-06-22 Thread Simon Sticko
Hi,

Lex, given your previous question:

https://groups.google.com/forum/#!topic/dealii/xghVE7Km78o

If you are already using an FEFaceValues object with a dummy-quadrature 
(created from FiniteElement::get_unit_face_support_points())
then you can use one of the following functions on FEFaceValues to get the 
location of the support points in real space:

quadrature_point()
get_quadrature_points()

Best,
Simon

On Monday, June 22, 2020 at 10:53:10 PM UTC+2, Lex Lee wrote:
>
> Hello Doug,
>
>
> Thanks for your help.
>
> However, I need to say, I just had done exactly what you described before 
> I posted this question.
>
> "FiniteElement::get_unit_face_support_points()” returns Points;
>
> “transform_unit_to_real_cell()” requires Point as input.
>
> The dimensions of points in cells is different from that on faces. That’s 
> one of the reasons why I posted this question.
>
>
> If I misunderstand you, please kindly let me know.
>
>
>
> Regards,
> Lex 
>
>
> On Jun 21, 2020, at 7:05 PM, Doug Shi-Dong  > wrote:
>
> Hello Lex,
>
> transform_unit_to_real_cell() takes a point as unit cell's point an input 
> and returns the physical point location. Therefore, you can give the unit 
> face support points from the FiniteElement::get_unit_face_support_points() 
> (assuming your FiniteElement has support points). And pass it to 
> transform_unit_to_real_cell().
>
> Doug
>
> On Sunday, June 21, 2020 at 8:10:30 PM UTC-4, Lex Lee wrote:
>>
>> Hello Doug,
>>
>> Thanks for your kind help.
>>
>> I am trying to understand your suggestion in the original email.
>>
>> Are you suggesting me using the member function 
>> "transform_unit_to_real_cell()" to get all support points in a unit cell? 
>> In the step, sort out the face points I need? Am I correct?
>>
>>
>> Regards,
>>
>> Lex
>>
>> On Sunday, June 21, 2020 at 5:00:35 PM UTC-7, Doug Shi-Dong wrote:
>>>
>>> Hello Lex,
>>>
>>> You were close.
>>>
>>>
>>> https://www.dealii.org/current/doxygen/deal.II/classMapping.html#ae5df63553eb8ed170c3b90524853dd48
>>>
>>> The project_real_point_to_unit_point_on_face() is useful is you take a 
>>> "volume" point within the cell. Since you are mapping from unit to real, 
>>> you would always know whether your point is on the face or not.
>>>
>>> Therefore, transform_unit_to_real_cell() should do the trick. Unless you 
>>> somehow actually need to project a "volume" point, but projecting within 
>>> unit cell is easy, which you would then follow with 
>>> transform_unit_to_real_cell().
>>>
>>> Best regards,
>>>
>>> Doug
>>>
>>> On Saturday, June 20, 2020 at 8:16:22 PM UTC-4, Lex Lee wrote:

 Hello Deal.II Users,


 I want to get the physical (geometry) coordinates for support points on 
 cell faces. Also I know that there are such member functions that can me 
 map points between reference and real cells in this link: 
 https://www.dealii.org/current/doxygen/deal.II/classMapping.html . 

 However, via the link mentioned above, I only find a member function 
 called project_real_point_to_unit_point_on_face 
 
  (real 
 to unit) , no functions like project_unit_point_to_real_point_on_face 
 (unit to real).  

 Could someone kindly tell me how I can solve this problem?  


 Regards,

 Lex

>>>
> -- 
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see 
> https://groups.google.com/d/forum/dealii?hl=en
> --- 
> You received this message because you are subscribed to a topic in the 
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit 
> https://groups.google.com/d/topic/dealii/uglzE4d4Flo/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to 
> dea...@googlegroups.com .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/fd8f078d-f0d4-4319-8977-a9a8e4672fb5o%40googlegroups.com
>  
> 
> .
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/aedf08b6-71cf-4452-9619-2155a497f35bo%40googlegroups.com.


[deal.II] Re: Error calculating shape gradient

2020-06-10 Thread Simon Sticko
Hi,
Yes. Here, you add a number of scalar elements of different degree:

> for (unsigned int deg=1; deg <= max_fe_degree; ++deg)
> {
> fe_collection.push_back(FE_Q(deg));
> quadrature_collection.push_back(QGauss(deg+1));
> face_quadrature_collection.push_back(QGauss(deg+1));
> }

which makes sense for a scalar problem, like the Poisson equation.

But if you solve for displacement in a material the solution has dim 
components:

u = (u_x, u_y, u_z)

so you want to use a FESystem element with dim components:

fe_collection.push_back(FESystem(FE_Q(deg), dim));

Best,
Simon


On Wednesday, June 10, 2020 at 12:03:57 PM UTC+2, A.Z Ihsan wrote:
>
> Hi Simon, 
>
> i followed the tutorial step-27 about hp-fem, using the similar 
> constructor with dim = 3
> for (unsigned int deg=1; deg <= max_fe_degree; ++deg)
>   {
> fe_collection.push_back(FE_Q(deg));
> quadrature_collection.push_back(QGauss(deg+1));
> face_quadrature_collection.push_back(QGauss(deg+1));
>   }
>
> I tried to look into the fe_collection.n_components(), it yields 1.
> is that what you meant?
>
> BR, 
> Ihsan
>
> On Wednesday, June 10, 2020 at 9:07:58 AM UTC+2, Simon Sticko wrote:
>>
>> Hi,
>> from the error message you can see that the element you are using only 
>> has 1 component. You get an error because you are trying to access 
>> component 1, which doesn't exist. Since your element should have dim 
>> components, there is likely something wrong with how you create your 
>> element. It should probably be similar as in step-8.
>>
>> Best,
>> Simon
>>
>>
>> On Wednesday, June 10, 2020 at 8:29:03 AM UTC+2, A.Z Ihsan wrote:
>>>
>>> Hi all, 
>>>
>>> i am getting this error while calculating the local cell matrix for my 
>>> hp fem application
>>>
>>> dealii::Tensor<1, spacedim> dealii::FEValuesBase>>>> spacedim>::shape_grad_component(unsigned int, unsigned int, unsigned int) 
>>>>> const [with int dim = 3; int spacedim = 3]
>>>>
>>>> The violated condition was: 
>>>>
>>>> component < fe->n_components()
>>>>
>>>> Additional information: 
>>>>
>>>> Index 1 is not in the half-open range [0,1).
>>>>
>>>> Stacktrace:
>>>>
>>>> ---
>>>>
>>>> #0  ../build/local: dealii::FEValuesBase<3, 
>>>>> 3>::shape_grad_component(unsigned int, unsigned int, unsigned int) const
>>>>
>>>> #1  ../build/local: dealii::SymmetricTensor<2, 3, double> 
>>>>> LinearElasticity::get_strain<3>(dealii::FEValues<3, 3> const&, unsigned 
>>>>> int, unsigned int)
>>>>
>>>> #2  ../build/local: LinearElasticity::HPSolver<3, 
>>>>> dealii::Vector 
>>>>> >::assemble_cell_matrix(dealii::TriaActiveIterator>>>> > 
>>>>> 3>, false> > const&, dealii::FullMatrix&, dealii::hp::FEValues<3, 
>>>>> 3>&)
>>>>
>>>> #3  ../build/local: 
>>>>> LinearElasticity::HPSerialSolver<3>::assemble_linear_system(LinearElasticity::HPSerialSolver<3>::LinearSystem&)
>>>>
>>>> #4  ../build/local: LinearElasticity::HPSerialSolver<3>::solve_problem()
>>>>
>>>>  
>>>>
>>>>
>>> The idea is similar to the tutorial step-27, but here i use the 
>>> symmetric tensor for returning strain. Can someone explain to me what the 
>>> error is and how to resolve it? 
>>>
>>> BR, 
>>> Ihsan
>>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/dac7f402-fb90-4f1c-940d-5450ddd0cf5do%40googlegroups.com.


[deal.II] Re: Error calculating shape gradient

2020-06-10 Thread Simon Sticko
Hi,
from the error message you can see that the element you are using only has 
1 component. You get an error because you are trying to access component 1, 
which doesn't exist. Since your element should have dim components, there 
is likely something wrong with how you create your element. It should 
probably be similar as in step-8.

Best,
Simon


On Wednesday, June 10, 2020 at 8:29:03 AM UTC+2, A.Z Ihsan wrote:
>
> Hi all, 
>
> i am getting this error while calculating the local cell matrix for my hp 
> fem application
>
> dealii::Tensor<1, spacedim> dealii::FEValuesBase>> spacedim>::shape_grad_component(unsigned int, unsigned int, unsigned int) 
>>> const [with int dim = 3; int spacedim = 3]
>>
>> The violated condition was: 
>>
>> component < fe->n_components()
>>
>> Additional information: 
>>
>> Index 1 is not in the half-open range [0,1).
>>
>> Stacktrace:
>>
>> ---
>>
>> #0  ../build/local: dealii::FEValuesBase<3, 
>>> 3>::shape_grad_component(unsigned int, unsigned int, unsigned int) const
>>
>> #1  ../build/local: dealii::SymmetricTensor<2, 3, double> 
>>> LinearElasticity::get_strain<3>(dealii::FEValues<3, 3> const&, unsigned 
>>> int, unsigned int)
>>
>> #2  ../build/local: LinearElasticity::HPSolver<3, dealii::Vector 
>>> >::assemble_cell_matrix(dealii::TriaActiveIterator>> > 
>>> 3>, false> > const&, dealii::FullMatrix&, dealii::hp::FEValues<3, 
>>> 3>&)
>>
>> #3  ../build/local: 
>>> LinearElasticity::HPSerialSolver<3>::assemble_linear_system(LinearElasticity::HPSerialSolver<3>::LinearSystem&)
>>
>> #4  ../build/local: LinearElasticity::HPSerialSolver<3>::solve_problem()
>>
>>  
>>
>>
> The idea is similar to the tutorial step-27, but here i use the symmetric 
> tensor for returning strain. Can someone explain to me what the error is 
> and how to resolve it? 
>
> BR, 
> Ihsan
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/9e6af2f7-695c-4a4f-bcd6-17c9a4f008cao%40googlegroups.com.


[deal.II] Re: How to get normal / tangential vectors at nodes, not at quadrature points?

2020-06-05 Thread Simon Sticko
Hi,
I suspect that the condition you have would be easier to enforce weakly by 
adding terms to the weak formulation, as is done in e.g. discontinuous 
Galerkin.

Nevertheless, if you still want the normals at the nodes you can create a 
"dummy" quadrature which has the nodes on the face as quadrature points:

const FiniteElement  = ... // some element 

const Quadrature quadrature(fe.get_unit_face_support_points());

You can then use this quadrature in FEFaceValues to get the normals. Here, 
the weights of the dummy quadrature are inf, but if the normals are 
everything you need that isn't a problem.

Best,
Simon

On Thursday, June 4, 2020 at 8:09:08 PM UTC+2, Lex Lee wrote:
>
> Hello Deal.II Users,
>
>
> I am working on setting a constraint: (V - v_s ) \cdot n =  \phi_f (v_f - 
> v_s) \cdot n  on three vector variables at the interface with affine 
> constraint class in Deal.II's library.
>
> (Where, V, v_s, and v_f are velocity vectors on two different domains,  
> and they are coupled with each other at the interface; n is the normal 
> vector; \phi_f is the volume fraction.)
>
>
> I decided to constrain either V_x= -(...) / n_x  or  V_y= -(...) / n_y 
> depending on the absolute values of n_x or n_y (n_x and n_y are the 
> components at x, y directions of the normal vector n).
>
> Now, I have a problem with getting the normal vectors *at nodes *when 
> playing with affine constraint within the fe_nothing and hp:finite element 
> framework. The relevant function/class I have found in Deal.II's library 
> all return the normal vector *at quadrature points,* not at nodes.
>
> Have you ever encountered a problem like this? Could you kindly share me 
> your idea of handling this problem? Thousand of thanks. 
>
>
> Best,
>
> Lei Li
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/ad68ea89-35dc-4f0d-a9a0-3eb2f2cedd59o%40googlegroups.com.


[deal.II] Re: Resolving with different RHS

2020-06-02 Thread Simon Sticko
Hi,
a better option than computing the inverse is to factorize the matrix. This 
can be done using the SparseDirectUMFPACK solver:

https://www.dealii.org/current/doxygen/deal.II/classSparseDirectUMFPACK.html

You might want to take a look at step 22, which uses this solver:

https://www.dealii.org/current/doxygen/deal.II/step_22.html

Best,
Simon

On Tuesday, June 2, 2020 at 12:04:38 PM UTC+2, Andreas Kyritsakis wrote:
>
> Dear all,
>
> I have a problem where a Poisson-like equation has to be solved again and 
> again in many time steps. At each timestep the LHS remains the same while 
> the RHS changes (slightly). My current implementation is to use the 
> standard SolverCG, passing the old solution as initial, which already 
> reduces the number of CG steps required. Yet, I wonder whether my approach 
> is naive and there is a faster way to implement this, taking better 
> advantage of the fact that the LHS stays always the same. I initially 
> thought about calculating the inverse matrix once and just doing matrix 
> multiplication, but since the mass matrix is a huge sparse matrix, only 
> storing the inverse (which in general is not sparse) would require huge 
> memory. I also thought about the LinearOperator concepts, but if I 
> understood correctly, they just implement a nice wrapper to call a solver 
> each time. Am I missing something?
>
> Cheers,
> Andreas
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/4bc6a09b-834c-485e-a26c-cbbbc6a45338%40googlegroups.com.


Re: [deal.II] Re: FE_Q

2019-11-30 Thread Simon Sticko
Hi,
Just to elaborate on what Daniel already said, Oliver, note that it is the 
constructor of the Step3-class that calls the constructor of FE_Q, 
that is

Step3::Step3()
  : fe(1)
  , dof_handler(triangulation)
{}

If you miss this call in the Step3-constructor, you will get the error you 
stated earlier.

Best,
Simon

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/633d2d99-c5fb-4720-94f4-38d4127d8b06%40googlegroups.com.


Re: [deal.II] Re: Bug Report in GridTools::get_active_neighbors()

2019-11-29 Thread Simon Sticko
Hi,
Fixed in #9108, already merged in. I should probably have written that here for 
completeness. Sorry.

Best,
Simon

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/eb765d3a-d9bc-4669-928d-01c281f73911%40googlegroups.com.


[deal.II] Re: Bug Report in GridTools::get_active_neighbors()

2019-11-28 Thread Simon Sticko
Hi,
I think it's just the documentation that can be slightly improved here.

The function get_active_child_cells has a similar signature but it has 
this note in its documentation:

* @note Since in C++ the MeshType template argument can not be deduced from
* a function call, you will have to specify it after the function name, as
* for example in
* @code
*   GridTools::get_active_child_cells > (cell)
* @endcode
*/
template 
std::vector
get_active_child_cells(const typename MeshType::cell_iterator );

The same note should maybe be added to get_active_neighbors for clarity.

Best,
Simon

On Thursday, November 28, 2019 at 10:03:17 AM UTC+1, Andreas Kyritsakis 
wrote:
>
> Hello,
>
> This is to report what I think is a bug in 
> GridTools::get_active_neighbors().
>
> The following code should compile but it does not.
>
> #include 
> #include 
> #include 
> #include 
> using namespace dealii;
>
> template
> class Test{
> public: 
> Test(Triangulation *tria) : triangulation(triangulation), 
> dof_handler(*tria){}
> 
> void run(){
> typename DoFHandler::active_cell_iterator cell;
>   
> for (cell = dof_handler.begin_active(); cell != dof_handler.end(); 
> ++cell){
> std::vector::active_cell_iterator> 
> neighbors;
> GridTools::get_active_neighbors(cell, neighbors); 
> }
> }
> private:
> Triangulation* triangulation;
> DoFHandler dof_handler;
> };
>
> int main(){
> Triangulation<2> triangulation;
> GridGenerator::hyper_cube(triangulation, -1, 1);
> Test<2> test();
> test.run();
> }
>
> I think somehow the compiler cannot resolve the generic type 
> MeshType::active_cell_iterator into the specific one 
> DoFHandler::active_cell_iterator, although from what I understand in 
> the documentation it should. Probably some template specification missing 
> somewhere, but I can't find it.
>
> In my code I fixed it easily by copying the function directly to my own 
> files and specifying appropriately, but I think it would be nice to report 
> the bug.
>
> Cheers,
> Andreas 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/0a78a60f-89d1-4c6a-9ba3-724d7112c355%40googlegroups.com.


Re: [deal.II] Re: fe_enriched and step-47

2017-04-09 Thread Simon Sticko
Hi.
Regarding this:

> On Thursday, April 6, 2017 at 1:51:41 PM UTC+2, Denis Davydov wrote:
> 
> The issue with integrating and quadrature would still remain.


I've been implementing an algorithm to generate bulk and surface quadrature 
given a level set function as boundary description. It's mainly an 
implementation of this algorithm: dx.doi.org/10.1137/140966290. 
I'm hoping to finish this within the coming weeks and contribute with this 
to dealii.

Best 
Simon

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.